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Preface 



On the occasion to submit my thesis to the preprint server, I would like to 
refer to previous and recent studies which are not referred to in the original 
thesis. 

In the early days of quantum annealing study Amara et al. used the 
imaginary time Schrodinger equation to search the global minimum state of 
the classical system [1]. Finnila et al. introduced a controlled and scheduled 
quantum effect to the Monte Carlo Simulation [2]. 

Recently, thermal and quantum annealing in the disordered Ising magnet, 
LiHoo.44Yo.56F4, with transverse magnetic field were compared[3,4]. Numer- 
ical studies of quantum annealing for the disordered Ising model and the 
protein folding problem were performed by a few different groups[5,6]. 

A quantum computer algorithm is also related with quantum annealing. 
The algorithm uses the adiabatic evolution of a time dependent Hamilto- 
nian[7,8,9,10]. 
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Abstract 



We introduce quantum fluctuations into the simulated annealing process of 
optimization problems, aiming at faster convergence to the optimal state. 
Quantum fluctuations cause transitions between states and thus play the 
same role as thermal fluctuations in the conventional approach. The idea 
is tested by the two models, the transverse Ising model and the traveling 
salesman problem. 

Adding the transverse field to the Ising model is a simple way to introduce 
quantum fluctuations. The strength of the transverse field is controlled as 
a function of time similarly to the temperature in the conventional method. 
The goal is to find the ground state of the diagonal part of the Hamiltonian 
with high accuracy as quickly as possible. We also consider the traveling 
salesman problem. This model can be described by the Ising spin, so that 
we also apply the same technique to the transverse Ising model. 

We solve the time-dependent Schrodinger equation numerically for small 
size systems with various types of exchange interactions of the Ising model. 
Comparison with the results of the corresponding classical (thermal) method 
reveals that the quantum method leads to the ground state with much larger 
probability in almost all cases if we use the same annealing schedule of the 
control parameters. We check the case of large-size systems by using the 
quantum Monte Carlo method. The simulation supports the results of the 
small-size systems, while the dynamics of the Schrodinger equation and the 
quantum Monte Carlo method are not the same. We find that the simu- 
lated annealing by quantum fluctuations has a better performance than the 
conventional method for the ground state search of the Ising-spin systems. 

The calculation of the traveling salesman problem is also performed as 
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an application to the general optimization problems. We obtain the same 
feature of the fast convergence to the optimal state as the transverse Ising 
model by using the quantum fluctuations. 

We also find that the relaxation time is quite short for quantum systems 
by numerical simulations. We consider this is one of the reasons why the an- 
nealing in quantum systems have a better performance of finding the optimal 
state in comparison with classical systems. 
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Chapter 1 
Introduction 



1.1 The Combinatorial Optimization Prob- 
lem 

The combinatorial optimization problem is to find a minimum or maximum 
value of a function of very many independent variables, where the variables 
take discrete values. This function, usually called the cost function or objec- 
tive function, represents a quantitative measure of the "goodness" of some 
complex systems. A famous example is the traveling salesman problem |], 0] 
to find the shortest or lowest-cost route to visit given cities. The techniques 
to find the shortest route have been studied for a long time, because such a 
problem actually happens in various situations, and people have to obtain the 
solution. One of the most important application of the optimization prob- 
lems today is the circuit design of the computer systems. The cost function 
is the wiring length in this case. Long wiring makes delays in electric circuit, 
and the speed of computer slows down. 

The algorithm to find the shortest tour of the traveling salesman problem 
or the wiring length of the circuit design in polynomial time of N, where N 
is the size of the problem, is not found yet. When N increases, the complex- 
ity of the calculation increases exponentially and goes over the limit of the 
computational power. This kind of problems are called the "Nondeterminis- 
tic Polynomial-time solvable hard or complete (iVP-hard or iVP-complete)" 
problems. (Details of the definitions are explained in the next section.) Any 
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algorithm of solving an iVP-hard or ./VP-complete problem in polynomial 
time is not found so far and people consider that there does not exist such 
an algorithm. Therefore we consider approximate algorithms by which we 
can obtain an approximate solution or an exact solution with some probabil- 
ity in polynomial time. In this thesis we focus on how to solve the iVP-hard 
or iVP-complete problems approximately as fast as possible. 



1.2 The iVP-hard and The iVP-complete Prob- 
lems 

In this section we explain several levels of difficulty of optimization problems. 
The iVP-hard and iVP-complete problems belong to the class quite hard to 
solve. The difference between the classes iVP-hard and iVP-complete is not 
in the degree of difficulty but the way to answer the questions; "yes or no" 
for iVP-complete and in general statements for iVP-hard. 

First, let us consider the class NP before the definition of the classes 
iVP-hard and ./VP-complete. The class NP is a subclass of the solvable class 
of decision problems whose answer is expressed by either "yes" or "no" . For 
instance, the Hamilton circuit problem is the problem to check the existence 
of the Hamilton circuit which is the closed loop constructed by edges and all 
the vertexes in a given graph. An example of the Hamilton circuit problem 
and its Hamilton circuit are shown in Fig. |1 . 1| . 

Roughly speaking, the NP problem is the problem which can be verified 
in polynomial time by the people who know the evidence to prove the exis- 
tence of the solution. For instance, one can answer "yes" for the Hamilton 



circuit problem shown in Fig. |1 . 1| , when the Hamilton circuit (the evidence 



in this case) is given as a bold line in Fig. |1 . 1| . The verification time of the 
evidence, checking the closed loop which contains all the vertexes, is pro- 
portional to N. To be precise, the NP problem is defined as a problem 
which can be solved by the nondeterministic Turing machine in polynomial 
time of N. However, the definition of the nondeterministic Turing machine 
is complicated, and the reader who is not interested in the detail can skip 
the following two paragraphs. 
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Figure 1.1: The Hamilton graph is the graph which has a Hamilton circuit. 
The bold line is the Hamilton circuit. 

Before presenting the definition of the nondeterministic Turing machine, 
we explain the definition of the usual (deterministic) Turing machine. The 
nondeterministic Turing machine is the extended model of the Turing ma- 
chine. The definition of the Turing machine is constructed with 

• the tape which includes alphabets (as commands and data). 

• the state of the machine. 

• the transition function of the state by using the alphabet read from the 
tape. 

The Turing machine is a simple model of actual computers. For the de- 
terministic Turing machine, the transition function is single-valued. Turing 
machine reads an alphabet from the tape and modifies its configuration (the 
state and alphabets on the tape) from the previous configuration according 
to the given alphabet. A unique modified state is given by this transition. 
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On the other hand, the transition function of the nondeterministic Tur- 
ing machine is multi-valued. The feature of the nondeterministic Turing 
machine quite differs from conventional computers. The configuration of the 
nondeterministic Turing machine becomes multiple configurations. This can 
be understood loosely in a way that the nondeterministic Turing machine 
makes replicas of possible configurations by a single step and the number 
of replicas increases exponentially as a function of the step. This feature 
is quite useful to solve complicated problems. For finding the ground-state 
energy of Ising-spin systems, for example, we have to enumerate all possible 
spin configurations, if the ground state is nontrivial. These configurations 
are illustrated by the tree structure in Fig. |1.2| . The direction of spins is 
assigned at each branching point and the level corresponds to the location of 
the spin. Each path means the configuration. The nondeterministic Turing 
machine can check all paths in parallel. At the first step, the machine makes 
a replica and the number of machines is two. Each machine makes a replica 
at each step so that the number of machines is 2 N at the iVth step. The 
calculations in the same level of the tree are performed at a time. By this 
procedure, we can enumerate 2 N configurations in iV steps. 

All NP problems can be verified in polynomial time even on a determin- 
istic Turing machine, but the difficulty of solving the problems is not the 
same. The problems are divided into some subclasses — the subclasses P for 
easy problems and ./VP-complete for difficult problems. The problem which 
has the algorithm to be solved in polynomial time by the deterministic Tur- 
ing machine belongs to the class of Polynomial-time solvable (P). Obviously 
NP contains P, because the nondeterministic Turing machine includes the 
deterministic Turing machine in its definition. 

For an example of the P problem, we explain the Hamilton circuit problem 
for the graph whose maximum degree of the edge in a vertex is limited to 
two belongs to P. The typical graphs are shown in Fig. |1.3| . The answers 
of the left and the right problems are "yes" and "no" respectively. One can 
check the existence of the Hamilton circuit by the following way: 

1. Start from a vertex and jump to the neighbor vertex connected by 
an edge. (One can not jump to the neighbor which has already been 




N spins 2 configurations 

Figure 1.2: The path means the spin configuration. The number of the 
leaves (the bottom ends of the tree) equals to the number of configurations 
as 2". 

visited once.) 

2. Repeat 1 until no vertex to jump is left. 

3. If the last vertex is connected to the starting vertex and all the vertexes 
have been visited, the answer is "yes" . The answer is "no" other wise. 

This procedure needs polynomial time. If the maximum degree of the edge 
in a vertex is not limited to two, this procedure does not work. 

Let us consider the difficult classes ./VP-hard and iVP-complete. If any 
problems in NP can be solved by using a particular algorithm of the problem 
repeatedly polynomial times at most, the problem is called "iVP-hard". In 
other wards, we can solve any iVP-problems by calling the subroutine for a 
particular iVP-hard problem (maybe once except for other iVP-hard prob- 
lems). The ./VP-hard problem is more difficult than any NP problems or as 
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Figure 1.3: The Hamilton circuit problems with the answer "yes" and "no". 
The left graph is the single loop (the Hamilton circuit) and right one is not 
the single loop. 

difficult as the other iVP-hard problems in polynomial order. The algorithm 
of the iVP-hard problem is the most general algorithm in NP class. The 
iVP-hard problems by definition are not limited to belong to the class NP. 
For instance, the traveling salesman problem is not an NP problem but an 
iVP-hard problem because the answer is not "yes or no" but the tour length 
and route. The subclass "./VP-complete" is the set of the ./VP-hard problems 
which belong to the class NP. The iVP-complete problems are more difficult 
than any other NP problems and as difficult as the iVP-hard problems. The 
Hamilton circuit problem is known as an ./VP-complete problem [J3J. 

The difficulty of the problems can be expressed symbolically as P < 
(NP - P - iVP-complete) < iVP-complete « iVP-hard. It is a trivial 
consequence of their definition that if one of the iVP-complete problems 
can be solved in polynomial time, then all the problems of the NP class 
are also polynomial. In this case one would have P = NP. So far no 
polynomial algorithm has been found for an iVP-complete problem, and the 
question of whether P is equal to NP is still open, although the general 
belief is that iVP-complete problems are not polynomial. The most probable 
situation is sketched in Fig. |1.4| . The NP problems can be divided into 
three subclasses, iVP-complete, P and the subclass which does not belong 
to the previous two subclasses. The class iVP-hard is separated into two 
subclasses, iVP-complete and the other part by whether the problem is the 
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decision problem or not. Empirically the time it takes to solve an ./VP-hard 
and an AP-complete problems tends to scale exponentially with the size N. 
All the problems can be solved by the enumeration method by which all 
the configurational space of the problem are enumerated and this calculation 
needs the time proportional to e N . 



NP 



NP-complete 



NP-hard 



Figure 1.4: The class of NP problems assuming that P ^ NP. Under this 
assumption it can be shown that there are NP problems that are neither 
iVP-complete nor P. 



1.3 Spin Glass Models and Their Relation to 
Optimization 

Combinatorial optimization problems appear in various disciplines. In sta- 
tistical mechanics, the energy is equivalent to the cost function of the op- 
timization problem. Until the 1970s, almost all the systems to be studied 
were homogeneous so that the energy landscape has a simple structure and 
the optimal solution (the ground state) is trivial. Although such a system 
has many variables, the degrees of freedom can be reduced to small numbers. 
For instance, spins turn to the same direction in the ground state of the fer- 
romagnetic model and the spin variables are reduced to single-spin variable. 
This reduction comes from the homogeneous interactions. In the 1970s the 
study of the spin glass started. The spin glass system is a non-homogeneous 
system and the ground state is nontrivial. The background of the physics 
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of spin glasses and its relation to the combinatorial optimization problem is 
explained in this section. 

The study of spin glasses started in 1972 when Cannella and Mydosh @] 
investigated the AC susceptibility of a dilute magnetic alloy Au-Fe and found 
out that the AC susceptibility has a nontrivial sharp cusp. In the alloy, the 
interaction between two spins of Fe atoms is expressed as the RKKY 
interaction || ||, |?J 

Jij oc -3- cos(2fc F ry ) , (1.1) 

where is the distance between two spins and is the Fermi wave number. 
From this interaction and a random distribution of Fe atoms in space, it is 
possible that the interaction is a positive or a negative random variable. 

Edwards and Anderson proposed the so-called Edwards-Anderson (EA) 
model || to introduce the effect of the RKKY interaction in a diluted mag- 
netic alloy. The EA model has nearest-neighbor interactions Jij. The value 
of the interaction distributes as 




where z is the number of nearest neighbors. The average of the distribution 
is Jq/z, and the standard deviation is Jj\fz. The Hamiltonian is given by 

H = —"^^JijSiSj . (1.3) 

n.n. 

The EA solution is not solved exactly. Sherrington and Kirkpatrick || inves- 
tigated an infinite-range model, the so-called SK model. The Hamiltonian is 
given by 

H = ~ Yl Ji i SiS i ~ h J2 S i> ( L4 ) 

<ij> i 

where < ij > denotes a summation over all spin pairs. This is the reason 
of the name of the infinite-range model. The summation is different from 
that in the EA model. If the average of distribution of interactions is equal 
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to zero and the external field is vanishing, a second order phase transition 
at a finite temperature T c {= J/k B ) occurs, and the spin glass phase realizes 
below T c . 

The statistical property as the SK model was almost solved by further 



studies [10, 11, O, 13, 14, 15, 16, 17, 18|, but it is difficult to analyze the EA 



model by the same approaches. The EA model is often studied by computer 
simulations, the Monte Carlo simulations or other methods. The great in- 
terest of the EA model or the ±J model whose interactions are distributed 
binary (at J and — J) is whether the same picture of the SK model — the 



picture of the ultrametricity and the structure of the P(q) [13|, 18] - - exists 
or not. One can check this picture from the calculation of the distribution 
function P(q) where the overlap q is calculated from the pure states which are 
the equilibrium states separated by high energy barriers and can be obtained 
by the Monte Carlo simulation. If the simulation starts from a random con- 
figuration which is equivalent to the configuration at high temperature, the 
system is stuck in a metastable state and it takes long time to escape from 
the metastable state because the basins of the pure states and the metastable 
states are separated by high energy barriers. To avoid this difficulty of the 
simulation, one can start the simulation from the ground state, which is close 
to the pure state at low temperature. However, finding the ground state is 
also difficult. The reason is the following: Spins on the some plaquettes in the 
EA model do not satisfy all the bonds illustrated in Fig. |1.5| . This situation 
occurs, when the product of all the bonds in the plaquette has a minus sign. 
The frustrated plaquettes are located at random so that we have to consider 
the global structure of the configuration which pays as small penalty of the 
energy as possible. Therefore the ground state is not the trivial configura- 
tion. The difficulties in studying the EA model are the slow relaxation and 
the non-trivial configuration of the ground state. The latter problem is a 
typical example of the combinatorial optimization problem. 

The SK model is connected with the bipartition problem directly. We 
divide the group into two parts. "Likes and dislikes" between the members 
is expressed as the numerical value. The task is to find the best partition of 
the members. The case of = 6 is illustrated in Fig. |1.6| . The members are 
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Figure 1.5: The frustrated plaquette. The solid and the zigzag lines are the 
ferromagnetic and antiferromagnetic bonds. 

divided into open and closed circles, and the black and the gray lines mean 
"likes" and "dislikes" between the members. Assigning +1 to the members 
in one part and —1 to the members in the other part, we can regard that the 
members are expressed by Ising spins {cr.;}. The cost function of the "likes 
and dislikes" between two members i and j can be obtained as Jij- 2 -^—- If 
the two members i and j are in the same part, "likes and dislikes" is counted 
as Jij. If the two members are in the different parts, the value is zero. We 
can calculate the total cost of "likes and dislikes" by summing up all the 
pairs as: 

Total Cost = J2 J v ai % +1 

(ii) 

= - Jij<Ji(?j + const. (1.5) 

(ij) 

This cost function is exactly the same as the Hamiltonian of the SK model. 
The bipartition problem and the finding the ground state of the SK model 
are known as the ./VP-hard problems [HJ [2U|] . 

1.4 Simulated Annealing 



As shown in the previous section, large computational time is needed to solve 
an iVP-complete or an iVP-hard problem exactly. The technique of simu- 
lated annealing (SA) was first proposed by Kirkpatrick et al. (^TJ as a general 
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Figure 1.6: The bipartition problems with N = Q. The open and closed 
circles are the divided members of the group. The black and the gray lines 
mean "likes" and "dislikes" between the members. 

method to solve optimization problems. This method is an approximate al- 
gorithm to find the optimal solution in finite time, but the solution converges 
to the optimal state in the infinite-time limit. The probability to find the 
ground state increase with time and converges to one in the infinite-time 
limit. For practical applications, it is important to find a sufficiently good 
solution, where "good" means that the cost function is close to the optimal 
value. 

The idea is to use thermal fluctuations to allow the system to escape 
from local minima of the cost function so that the system reaches the global 
minimum under an appropriate annealing schedule (the rate of decrease of 
temperature). Let us consider a simple example of a system described by 
Ising spins {<7j ± 1}, binary variables. Spins interact each other as —JoiO^ 
and the total energy (cost function) of this system is H = — J a i a ji where 
ij runs over all possible pairs of sites which interact with each other. If the 
two spins <7j, <jj take the same value, the energy becomes lower. This system 
is the simplest model of a ferromagnet. In the optimal configuration, all spins 
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align in the same direction, up or down. 

Once we define the energy of the system, we can consider statistical me- 
chanics of this system. The spin configuration {<Ji} is modified by thermal 
noise. The probability with which the configuration {a{\ appears is propor- 
tional to the Boltzmann factor exp(—E({<Ji}) / k-gT) , where T is the tempera- 
ture and k-Q is the Boltzmann constant. The probability has to be normalized 
and the normalization factor is 1/Z. The inverse of the factor, the so-called 
partition function, is expressed as: 

Z= exp(-£({a,})/A;BT) (1.6) 

{<Ti=±l} 

The lower the temperature becomes, the larger the probability of the 
spin configuration of the lowest energy becomes. The system freezes to the 
ground state as the temperature goes to zero. However, the system does 
not always freeze to the ground state in rapid cooling. As illustrated in 
Fig. |1.7| , the system may be trapped in a local minimum. Therefore, the 
temperature should be cooled attentively. At high temperature, the system 
changes the configuration almost randomly and searches the global structure 
of the energy landscape. The system seeks the area in which the valley of 
the global minimum is included when the temperature is decreased. Finally, 
the system is stuck at the valley of the global minimum at zero or sufficiently 
low temperature (lower than the energy gap of the ground state and the first 



excited state). The idea of this cooling is illustrated in Fig. |L8 

By the way, how can we cool the system slowly enough? Geman and 
Geman proved a theorem on the annealing schedule for combinatorial op- 



timization problems piq| . They showed that any system reaches the global 
minimum of the cost function asymptotically if the temperature is decreased 
as T = c/logt or slower, where c is a constant determined by the system 
size and other structures of the cost function. When the product of the tem- 
perature and the Boltzmann constant k-^T become less than the energy gap 
between the ground state and the first excited state, the system jumps to 
the excited state rarely by thermal fluctuation. Writing such a temperature 
as Tq, we obtain the time to find the ground state as t ~ exp(c/ k-gTo) . This 
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Configuration 

Figure 1.7: Cooling from high temperature to low temperature. If the 
cooling rate is too fast, the system may be trapped in a local minimum. 

bound on the annealing schedule may be the optimal one under generic con- 
ditions although faster decrease of the temperature often gives satisfactory 
results in practical applications for many systems. 

1.5 Quantum Annealing 

In SA thermal fluctuations are introduced in the optimization problem so 
that transitions between states take place in the process of search for the 
global minimum among many states. Thus there seems to be no reasons 
to avoid use of other mechanisms for state transitions if these mechanisms 
may lead to better convergence properties. One such possibility is the gen- 
eralized transition probability of Tsallis lf23|| , which is a generalization of the 
conventional Boltzmann-type transition probability appearing in the master 
equation and thus used in Monte Carlo simulations. In this method, the 
system converges to the optimal state under power-low decrease of the tern- 
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Configuration 

Figure 1.8: The system seek the various levels of the energy-landscape 
structure corresponding to the temperature. The system searches the global 
(local) structure at high (low) temperatures. 



perature |23|, [24|, gg, |26[] . This schedule is faster than the log schedule of the 



conventional SA. However, the generalized transition probability does not 
satisfy the condition of the detailed balance, the sufficient condition for con- 
vergence to the equilibrium state, and the physical meanings of the system 
is not trivial. 

We seek another possibility of making use of quantum tunneling processes 
for state transitions, which we call quantum annealing (QA). In particular 
we would like to learn how effectively quantum tunneling processes possibly 
lead to the global minimum in comparison to temperature-driven processes 
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used in the conventional method of SA. 

In the virtual absence of previous studies along such a line of consid- 
eration, it seems better to focus our attention on a specific model system, 
rather than to develop a general argument, to gain an insight into the role 
of quantum fluctuations in the situation of optimization problem. Quantum 
effects have been found to play a very similar role to thermal fluctuations in 
the Hopfield model and the image restoration problem in a transverse field 



in thermal equilibrium |27|, g^j . This observation motivates us to investigate 
dynamical properties of the Ising model under quantum fluctuations in the 
form of a transverse field. We therefore start the study of the QA by the 
transverse Ising model with a variety of exchange interactions. The trans- 
verse field controls the rate of transition between states and thus plays the 
same role as the temperature does in SA. We assume that the system has no 
thermal fluctuations in the QA context and the term "ground state" refers 
to the lowest-energy state of the Hamiltonian without the transverse field 
term. 

Static properties of the transverse Ising model have been investigated 
quite extensively for many years |3U], pHf . There have however been very 
few studies on the dynamical behavior of the Ising model with a transverse 
field. We refer to the work by Sato et al. who carried out quantum Monte 
Carlo simulations of the two dimensional EA model in an infinitesimal trans- 
verse field, showing a reasonably fast approach to the ground state |32|. We 
present a new point of view by comparing the efficiency of QA directly with 
that of classical SA in reaching the ground state. 

In the next chapter, we consider the time-dependent Schrodinger equa- 
tion for QA and the classical master equation for SA numerically for small- 
size systems with the same exchange interactions under the same annealing 
schedules. We explain the model and define the measure of closeness of the 
system in QA to the desired ground state. This measure is compared with 
the corresponding classical probability that the system is in the ground state. 
Calculations of probabilities that the system is in the ground state at each 
time for both classical and quantum cases give important implications on 
the relative efficiency of the two approaches. The numerical results for QA 
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and SA suggest that QA generally gives a larger probability to find to the 
ground state than SA under the same conditions on the annealing schedule 
and interactions. We also calculate the analytical solutions for the one-spin 
case which turns out to be quite non-trivial. Explicit solutions yield very 
useful information to clarify several subtle aspects of the problem. 

In Chapter 3, we extend our method by using the Monte Carlo method 
to calculate large-size problems. We adopt the master equation as the time- 
evolution equation of SA in contrast to the Schrodinger equation of QA. The 
original definition of SA is expressed in terms of the Monte Carlo method, 
which is equivalent to the master equation. It is natural to consider the 
applicability of Monte Carlo method to QA. The time-evolution of the wave 
function is expressed as: 

|tf(t)) = Te-^> ww |^(0)), (1.7) 

where the symbol T is the time ordering operator. There are two approaches 
to take into account the quantum effects or the time evolution gradually. The 
first one is to divide the Hamiltonian into some parts which are easy to diag- 
onalize (for instance, into three terms which include a x , a y , a z , respectively) 

as 

e n = lim ( e Ho,m e Hl/m •••)"• (1-8) 

This division is the Suzuki- Trotter decomposition. By inserting the complete 
set ^ \i)(i\(= 1), the system is mapped to a (d + l)-dimensional classical 
system. Once the system is expressed as a classical system, we can perform 
the Monte Carlo method, the so-called quantum Monte Carlo method. The 
second one is the division of the short-time evolution as 

e -i^dtH(t) = j im e -*AtW(t) . . . e -*AtW(0) /j gx 

At^O ' 

the path integral Monte Carlo method. The dynamics of the two approaches 
are not exactly the same as the dynamics of the Schrodinger equation. The 
performances of the two approaches are the same or better in comparison with 
the dynamics of the Schrodinger equation in some sense (a detailed discussion 
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is shown later). We finally adopt the quantum Monte Carlo method for the 
calculation of QA, because this method is similar to the original SA and we 
can find the optimal state precisely by the time evolution. The results of 
the simulations show that QA performs in finding the optimal state better 
than SA in spite of the disadvantage in the calculation of the quantum effect 
(the calculation on (d+ l)-dimensional systems is larger than (^-dimensional 
system). We also calculate the difference between the quenched system and 
annealed system. The calculation of the quenched system is almost the same 



as the study by Sato et al. [02| . In general, quenched systems tend to be 
stuck in local minima. We find that the classical quenched system is stuck 
more often than the quantum quenched system. We consider this is one of 
the reasons why QA works better than SA. 

So far, the systems are physical models. The problems of the ground 
state finding of physical models occupy the small area of the combinatorial 
problems and we have to check applicability to the general problems. One of 
the well-known problems is the TSP. We apply QA to the traveling salesman 
problem (TSP) in Chapter 4. TSP can be mapped to the Ising spin system 
so that we can use the same framework of QA. The modification of the 
configuration of TSP in SA is not the same as the modification of the usual 
one-spin dynamics of SA. Four spins are changed at a time in TSP. The 
quantum-effect term of the Hamiltonian (the cost function) has to be changed 
to the product of the four spin operators to fit the dynamics of SA. However, 
we adopt the same quantum-effect term, the transverse field, because of the 
difficulty of the calculation. In spite of this, QA also works and improves 
in comparison to SA. This fact suggests that QA can be applied to various 
problems in which we can not define the quantum effect exactly. 



Chapter 2 

Analysis of Small-size Systems 
by Differential Equations 

In this chapter, we introduce the quantum annealing (QA) for small-size sys- 
tems by using the Schrodinger equation and compare the results of QA and 
the simulated annealing (SA). The models (problems) we deal with here are 
the models in statistical physics, precisely the Ising model. The transverse 
field is introduced as the quantum effect and scheduled as a decreasing func- 
tion of the time. This field flips single spin at a time so that the transition 
is quite similar to the usual single-spin-flip dynamics of SA. 

The analysis is performed numerically and analytically for SA and QA. 
The system size of the numerical calculation is up to eight. The calculation 
limit of the size is around fourteen by the middle-class supercomputer today, 
(for a fourteen-spin system, calculation needs the storage about 2G Byte!) 
The models we use as test-bed are ferromagnetic, frustrated and random- 
interaction models. We also study for the single-spin quantum system with 
longitudinal and transverse field. The exact solutions are obtained for some 
schedules. The single-spin system is too small to discuss the property of QA, 
but the result provides some information for the results of the numerical 
calculations. 
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2.1 The Transverse Ising Model 

Let us consider the following Ising model with longitudinal and transverse 
fields: 

nit) = -^J^J-/iJ>f-r(*)5>f (2.i) 

ij i i 

ee ft -r(0j>*, (2.2) 

i 

where the type of interactions will be specified later. The term of longitu- 
dinal field was introduced to remove the trivial degeneracy in the exchange 
interaction term coming from the overall up-down symmetry that effectively 
reduces the available phase space by half. The T(t)-term causes quantum 
tunneling between various classical states (the eigenstates of the classical 
part 'Ho)- By decreasing the amplitude T(t) of the transverse field from a 
very large value to zero, we hopefully drive the system into the optimal state, 
the ground state of Tio. 

The natural dynamics of the present system is provided by the Schrodinger 
equation 

.« = «WI«. (2.3) 

We solve this time-dependent Schrodinger equation numerically for small-size 
systems. The representation to diagonalize Ho (the ^-representation) will be 
used hereafter. 

The corresponding classical SA process is described by the master equa- 
tion 

^ = $>,P,(t), (2.4) 

3 

where Pi(t) represents the probability that the system is in the ith state. We 
consider single-spin flip processes with the transition matrix elements given 

as 

exp(-Ei/T(t)) 



exp(-Ei/T(t)) + eM-Ej/T(t)) 
(otherwise) 



single-spin difference) 



-J2k^ c ki (i = j) ' ( 2 ' 5 ^ 



CHAPTER 2. ANALYSIS BY DIFFERENTIAL EQUATIONS 20 



In SA, the temperature T(t) is first set to a very large value and then 
is gradually decreased to zero. The corresponding process in QA should be 
to change T(t) from a very large value to zero. The reason is that the high- 
temperature state in SA is a mixture of all possible states with almost equal 
probabilities, and the corresponding state in QA is the linear combination of 
all states with equal amplitude in the ^-representation, which is the lowest 
eigenstate of the Hamiltonian Q2.1|) for very large T. The low-temperature 
state after a successful SA is the ground state of Ho, which should also be the 
eigenstate of H(t) as T(t) is reduced to zero sufficiently slowly in QA. Another 
justification of identification of T and T comes from the fact that the T = 
phase diagram of the Hopfield model in a transverse field has almost the 
same structure as the equilibrium phase diagram of the conventional Hopfield 
model at finite temperature if we identify the temperature axis of the latter 



phase diagram with the T axis in the former [27]. We therefore change T(t) in 
QA and T(t) in SA from infinity to zero with the same functional forms T(t) = 
T{t) = c/t, c/yt, cj log(t + 1) (t : — ► oo) or — ct (t : — oo — > 0). The reason 
for choosing these functional forms are that either they allow for analytical 
solutions in the single-spin case as shown in Sec. |2.3| or for comparison with 
the Geman-Geman bound mentioned in the previous chapter. 

To compare the performance of the two methods QA and SA, we calculate 
the probabilities P QA (t) = \(g\^(t))\ 2 for QA and P SA (t) = P g (t) for SA, 
where P g (t) is the probability to find the system in the ground state at time 
t in SA and \g) is the ground-state wave function of 7io- Note that we treat 
only small-size systems (the number of spins N = 8) and thus the ground 
state can be picked out explicitly. In the ideal situation Pqa(£) and Psa(£) 
will be very small initially and increase towards 1 as t — > oo. 

It is useful to introduce another set of quantities P| A (T) and Pq A (T). The 
former is the Boltzmann factor of the ground state of Tio at temperature T 
while the latter is defined as K^l^r)] 2 , where the wave function ip^ is the 
lowest-energy stationary state of the full Hamiltonian ( [2.1] ) for a given fixed 
value of T. In the quasi-static limit, the system follows equilibrium in SA 
and thus P SA (t) = P§l{T{t)). Correspondingly for QA, P QA (t) = P| A (r(t)) 
when r(t) changes sufficiently slowly. Thus the differences between both 
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sides of these two equations give measures how closely the system follows 
quasi-static states during dynamical process of annealing. 

2.2 Numerical results 

We now present numerical results on Psa and Pqa for various types of ex- 
change interactions and transverse fields. All calculations were performed 
with a constant longitudinal field h = 0.1 to remove trivial degeneracy. 

2.2.1 Ferromagnetic Model 

Let us first discuss the ferromagnetic Ising model with J = const for all pairs 



of spins. Figure pTT] shows the overlaps for the case of T(t) = T{t) = 3/ log(i+ 
1). It is seen that both QA and SA follow stationary (equilibrium) states 
during dynamical processes rather accurately. In SA the theorem of Geman 
and Geman [Q] guarantees that the annealing schedule T(t) = c/ log(l + 1) 
assures convergence to the ground state (Psa — > 1 in our notation) if c 
is adjusted appropriately. Our choice c = 3 is somewhat arbitrary but the 
tendency is clear for Psa — > 1 as t — > oo, which is also clear from approximate 
satisfaction of the quasi-equilibrium condition Psa(^) = PsA^Tif)). Although 
there are no mathematically rigorous arguments for QA corresponding to 
the Geman-Geman bound, the numerical data indicate convergence to the 
ground state under the annealing schedule Y{t) = 3/log(t + 1) at least for 
the ferromagnetic system. It should be remembered that the unit of time 
is arbitrary since we have set H — 1 in the Schrodinger equation (|2.3j) and 



the unit of time r = 1 in the master equation (|2.4j ). Thus the fact that the 



curves for QA in Fig. |2.1| lie below those for SA at any given time does not 
have any positive significance. 

If we decrease the transverse field and the temperature faster, T(t) = 
T(t) = 3/y/i, there appears a qualitative difference between QA and SA as 



shown in Fig. [2.2| . The quantum method clearly gives better convergence to 
the ground state while the classical counterpart gets stuck in a local minimum 
with a non-negligible probability. To see the rate of approach of Pqa to 1, 
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Figure 2.1: Time dependence of the overlaps Psa(^), -Pqa(^), -Psa(^W) an d 
P§ A (T(t)) of the ferromagnetic model with T(t) = T(t) = 3/log(t + 1). 



we have plotted 1 — Pqa in a log-log scale in Fig. pT3[ It is seen that 1 — Pqa 
behaves as const/t in the time region between 100 and 1000. 

By a still faster annealing schedule T(t) = T(t) = 3/t, the system becomes 



trapped in intermediate states both in QA and SA as seen in Fig. EOl. 



2.2.2 Frustrated Model 



We next analyze the interesting case of a frustrated system shown in Fig. [275 . 
The full lines indicate ferromagnetic interactions while the broken line is for 
an antiferromagnetic interaction with the same absolute value as the ferro- 
magnetic ones. If the temperature is very high in the classical case, the spins 
4 and 5 are changing their states very rapidly and hence the effective inter- 
action between spins 3 and 6 via spins 4 and 5 will be negligibly small. Thus 
the direct antiferromagnetic interactions between spins 3 and 6 is expected to 
dominate the correlation of these spins, which is clearly observed in Fig. ^15] 
as the negative value of the thermodynamic correlation function (o"|crg) c in 
the high-temperature side. At low temperatures, on the other hand, the spins 
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Figure 2.2: Time dependence of the overlaps of the ferromagnetic model with 

r(t) = T{t) = 3/Vt. 
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Figure 2.3: Time dependence of 1 — Pqa(£) of the ferromagnetic model with 
r(i) = 3/y/i. The dotted line represents t^ 1 to guide the eye. 



CHAPTER 2. ANALYSIS BY DIFFERENTIAL EQUATIONS 24 




Figure 2.4: Time dependence of the overlaps of the ferromagnetic model with 
r(t) = T{t) = 3/t. 

4 and 5 tend to be fixed in some definite direction and consequently the ef- 
fective ferromagnetic interactions between spins 3 and 6 are roughly twice 
as large as the direct antiferromagnetic interaction. This argument is justi- 
fied by the positive value of the correlation function at low temperatures in 
Fig. |2.6| . Therefore the spins 3 and 6 must change their relative orientation 
at some intermediate temperature. This means that the free-energy land- 
scape goes under significant restructuring as the temperature is decreased 
and therefore the annealing process should be performed with sufficient care. 

If the transverse field in QA plays a similar role to the temperature in 
SA, we expect similar dependence of the correlation function (cferf} g on the 
transverse field T. Here the expectation value is evaluated by the stationary 
eigenfunction of the full Hamiltonian ( EO ) with the lowest eigenvalue at 
a given T. The broken curve in Fig. |2.6| clearly supports this idea. We 
therefore expect that the frustrated system of Fig. |2.5| is a good test ground 
for comparison of QA and SA in the situation with a significant change of 
spin configurations in the dynamical process of annealing. 
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Figure 2.5: The frustrated model where the solid lines denote ferromagnetic 
interactions and the broken line is for an antiferromagnetic interaction. 
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Figure 2.6: Correlation functions of spins 3 and 6 in Fig. ^75] for the classical 
and quantum cases. In the classical model (full line) the correlation is shown 
as a function of temperature while the quantum case (dotted line) is regarded 
as a function of the transverse field. 
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The results are shown in Fig. &7\ for the annealing schedule T(t) = T(t) = 
3/^/i. The time scale r is normalized as r = £T C 2 in SA and r = tT^ in QA. 
The values, T c and T c , are the points where the correlation functions vanish 



in Fig. |2.6| . Thus both classical and quantum correlation functions vanish at 
t = 1. The tendency is clear that QA is better suited for ground-state search 
in the present system. 



0.6 



. 4 



0.2 




_j i i 



10 



100 



1000 



Figure 2.7: Time dependence of the overlaps of the frustrated model under 
T(t) = T(t) = 3/y/t. Here the time scale r is normalized by T c and T c (the 
points where the correlation functions vanish in Fig. |2.6| .) 



2.2.3 Random Interaction Model 

The third and final example is the Sherrington-Kirkpatrick (SK) model of 
spin glasses ||. Interactions exist between all pairs of spins and are chosen 
from a Gaussian distribution with vanishing mean and variance 1/N (N = 8 



in our case). Figure |2.8| shows a typical result on the time evolution of 
the probabilities under the annealing schedule T(t) = T(t) = 3/y/i. We 
have checked several realizations of exchange interactions under the same 
distribution function and have found that the results are qualitatively the 
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same. Figure again suggests that QA is better suited than SA for the 
present optimization problem. 
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Figure 2.8: Time dependence of the overlaps for the SK model with T(t) = 
T(t) = 3/Vi. 



2.3 Solution of the single-spin problem 

It is possible to solve the time-dependent Schrodinger equation explicitly 
when the problem involves only a single spin and the functional form of the 
transverse field is T(t) = —ct,c/t or c/y/i. We note that the single-spin 
problem is trivial in SA because there are only two states involved (up and 
down) and thus there are no local minima. This does not mean that the 
same single-spin problem is also trivial in the quantum mechanical version. 
In QA with a single spin, the transition between the two states is caused 
by a finite transverse field. The system goes through tunneling processes to 
reach the other state, and an appropriate annealing schedule is essential to 
reach the ground state. On the other hand, in SA, the transition from the 
higher state to the lower state takes place even at T = and thus the system 
always reaches the ground state. 
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Let us first discuss the case of T(t) = —ct with t changing from — oo to 
0. This is the well-known Landau-Zener model and the explicit solution of 
the time-dependent Schrodinger equation is available in the literature 



34], H, H, H. With the notation a(t) = (+\ip(t)) and b(t) = (-\ip(t)) and 
the initial condition o(— oo) = b(—oo) = l/y/2 (the lowest eigenstate), the 
solution for b(t) is found to be (see Appendix) 

, . , he- nh2 / 8c f 2ct + h , . ih 2 + 2c 

(2.6) 

where -D_a-i ; -D-a-2 represent the parabolic cylinder function (or Weber 
function) and z and A are given as 



2ce- 7Ti/ H, (2.7) 
ih 2 

= To- (2 ' 8 > 

The final value of b(t) at t = is 

= _twnz e 

2^ 

J 1 ^-(1 + ^/20) 1 

\r(l + i/i 2 /4c) hr(3/2 + z/i 2 /4c) J ' 

The probability to find the system in the ground state at t = is, when 
fc 2 /c» 1, 

^Qa(O) = |a(0)| 2 = 1 - |6(0)| 2 -1-^. (2-10) 

Thus the probability Pqa(^) does not approach 1 for finite c. 

We next present the solution for T(t) = c/t with t changing from to oo 
under the initial condition a = b = l/y/2 (see Appendix): 

b(t) = - 7 =e <w * fc F(l + ic, 1 + 2ic; -2iht), (2.11) 
v2 

where F is the confluent hypergeometric function. The asymptotic form of 
b(t) as t — > oo is 

^ ^(2^-^(2^) x ^^-^72 + ce ^+W2 (2M) -i} ( 2.i2) 
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The probability to find the system in the target ground state behaves asymp- 
totically as 

P QA (t) = \a(t)\ 2 (2.13) 

= l-l&WI 2 (2-14) 
sinh(7rc) 



1 - 



sinh(27rc) 



f wr ccos(2ht) cV c ] , rt . 

~ l-e- 2nc , (2.16) 

the last approximation being valid for c ^> 1 after t — > oo. The system does 
not reach the ground state as t — > oo as long as c is finite. Larger c gives 
more accurate approach to the ground state, which is reasonable because it 
takes longer time to reach a given value of T(= c/t) for larger c, implying 
slower annealing. 

The final example of the solvable model concerns the annealing schedule 
T(t) = c/\/i. The solution for b(t) is derived in Appendix under the initial 
condition a = b = 1/ \f2 as 



b{t) = ^e M F^-i>y,l-2iht) 



where 7 = c 2 /2h. The large-t behavior is found to be 



b(t) ~ y/ne 



-7TC 2 /4h 



1 A /77p(5/4)7Tj 

^V2r(|-*7) cT(-<7) 



f e -(V4)« 

+e i/li (2/ i t)^ 1 / 2+ ^ { — + 



.2 

c 



\ V2r(^ 7 ) 2v^r(i + 27) 

and the probability Pqa(oo) for c 2 /h 3> 1 is obtained as 



(2.18) 



P QA (oc) = l-|6(oo)| 2 ~l-^. (2.19) 
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This equation indicates that the single-spin system does not reach the ground 
state under the present annealing schedule T(t) — c/y/t for which the nu- 
merical data in the previous section suggested an accurate approach. We 
therefore conclude that the asymptotic value of -Pqa(0 in the previous sec- 
tion may not be exactly equal to 1 for T(t) — S/y/i although it is very close 
to 1. 

The annealing schedule T(t) — c/y/t has a feature which distinguishes 
this function from the other ones — ct and c/t. As we saw in the previous 
discussion, the final asymptotic value of -PqaCO is not 1 if the initial condition 
corresponds to the ground state for T — > oo, a = b = l/y/2. However, as 
shown in Appendix, by an appropriate choice of the initial condition, it is 
possible to drive the system to the ground state if T(t) = c/y/t. This is not 
possible for any initial conditions in the case of T(t) = —ct or c/t. 

2.4 Summary 

In this chapter we have started the study of quantum annealing (QA) from 
the transverse Ising model obeying the time-dependent Schrodinger equa- 
tion. The transverse field term was controlled so that the system approaches 
the ground state. The numerical results on the probability to find the sys- 
tem in the ground state were compared with the corresponding probability 
derived from the numerical solution of the master equation representing the 
SA processes. We have found that QA shows convergence to the optimal 
(ground) state with larger probability than SA in all cases if the same an- 
nealing schedule is used. The system approaches the ground state rather 
accurately in QA for the annealing schedule Y = c/y/t but not for faster 
decrease of the transverse field. 

We have also solved the single-spin model exactly for QA in the cases of 
Y(t) = —ct,c/t and c/y/t. The results showed that the ground state is not 
reached perfectly for all these annealing schedules. Therefore the asymptotic 
values of -Pqa(^) in numerical calculations are probably not exactly 1 although 
they seem to be quite close to the optimal value 1. 

The rate of approach to the asymptotic value close to 1, 1 — Pqa(£), 
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was found to be proportional to 1/t in Fig. |2.3| for the ferromagnetic model. 
On the other hand, the single-spin solution shows the existence of a term 



proportional to see Eq. (|2.18|) . Probably the coefficient of the 1/yt- 

term is very small in the situation of Fig. |2.3| and the next-order contribution 
dominates in the time region shown in Fig. |2.3| . 

A simple argument using perturbation theory yields useful information 
about the asymptotic form of the probability function if we assume that the 
system follows quasi-static states during dynamical processes. The probabil- 
ity to find the system in the ground state is expressed using the perturbation 
in terms of 1) as 

PQA(r)~i-r"£ 1 , (2.20) 

AO) 



where E\ is the energy of the ith state of the non-perturbed (classical) 
system and Eq is the ground-state energy. If we set T = c/ we have 

^-^sUv)'- (2 - 2i) 

Thus the approach to the asymptotic value is proportional to 1/t as long as 
the system stays in quasi-static states. The corresponding probability for SA 
is 

p -Eo/T 

^sa(T) ~ ^— ^ ~ 1 - X)e-^"^, (2.22) 

which shows absence of universal (1/t-like) dependence on time. 

The present method of QA bears some similarity to the approach by the 
generalized transition probability in which the dynamics is described by the 
master equation but the transition probability has power-law dependence on 
the temperature in contrast to the usual exponential form of the Boltzmann 
factor [ 23" . This power- law dependence on the temperature allows the system 
to search for a wider region in the phase space because of larger probabilities 
of transition to higher -energy stcites £tt a given T(t), which may be the reason 
of faster convergence to the optimal states [^3], [26| . The transverse field term 
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T in our QA represents the rate of transition between states which is larger 
than the transition rate in S A (see ( |2.5| ) ) at a given small value of the control 
parameter T(t) = T(t). This larger transition probability may lead to a 
more active search in wider regions of the phase space, leading to better 
convergence similarly to the case of the generalized transition probability. 

We have solved the Schrodinger equation and the master equation di- 
rectly by numerical methods for the purpose of comparison of QA and SA. 
This method faces difficulties for larger N because the number of states 
increases exponentially as 2^. The classical SA solves this problem by ex- 
ploiting stochastic processes, Monte Carlo simulations, which have the com- 
putational complexity growing as a power of N. The corresponding reduction 
of the computational complexity is lacking in QA, because the Schrodinger 
equation is not replaced directly by the stochastic processes. While the 
quantum Monte Carlo is not equivalent to the Schrodinger equation, we will 
conduct QA in the quantum Monte Carlo framework in the next chapter. 
Another problem is to devise implementations of QA in other optimization 
problems such as the traveling salesman problem or the Hamilton problem. 
The implementation of QA to the traveling salesman problem is explained 
as an example in Chapter f|. 



Chapter 3 



Monte Carlo Analysis of Larger 
Systems 

We will apply the Monte Carlo method to QA in this chapter. Two possible 
ways to solve QA are considered. Comparison between SA and QA by the 
Monte Carlo method is also performed. 

We desire to solve SA and QA by another method, because the calcula- 
tions of the master equation for SA and the Schrodinger equation for QA be- 
come difficult for large-size systems. As shown in the previous chapter, when 
we search the ground state by SA and QA, the operation of 2 7V -dimensional 
matrices is needed to solve the differential equations in the Ising-spin sys- 
tems. For example, real and complex 1024 x 1024 matrixes are needed for SA 
and QA with ten spins, respectively. To avoid such a difficulty of solving the 
master equation, the calculation is performed by the Monte Carlo method in 
the original definition of SA. Solving SA by the master equation is the equiv- 
alent to the Monte Carlo method. From the relation of the master equation 
and the Monte Carlo method in SA, we consider QA by the Monte Carlo 
method. 

Another motivation of considering the Monte Carlo method is the fol- 
lowing: In the actual procedure of QA, all the elements of the Hamiltonian 
matrix are calculated to solve the the Schrodinger equation. If the matrix is 
expressed in the ^-representation, the diagonal elements are the classical-term 
energy of each configuration. This implies that there is no advantage to the 
enumeration method (to enumerate all the possible configurations), because 
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we already calculated all the energy of possible configurations. Therefore, we 
have to solve QA by another method which requires less calculations than to 
solve the differential equation, like the Monte Carlo method. 

We consider the two methods, the path-integral Monte Carlo and the 
quantum Monte Carlo methods, for QA. The results of the two methods are 
not the same as the result of the Schrodinger equation, but the ground state 
can be found by both of the methods. We adopt the quantum Monte Carlo 
method for QA and compare the results of SA and QA by the Monte Carlo 
simulation. 



3.1 Monte Carlo Method 

The Monte Carlo method is powerful to analyze various problems of statis- 
tical mechanics. In classical statistical mechanics, we obtain the expectation 
value of a quantity A as a function of the temperature T = 1/(3 in the 
following form: 

where the summation runs over all the possible configurations of the system 
and E corresponds to the energy of each configuration. The number of all 
possible configurations is S , if the system includes N sites and each site 
takes S states independently. The number of the terms in the summation 
increases exponentially as a function of the system size N. It is difficult to 
calculate the whole summation for large-size systems. 

We consider replacing the summation with the weighted sampling. Equa- 
tion (R.ll) is rewritten as 



where Pi = e~@ El / e~^ E i . Pj is the Gibbs distribution, the probability of 
the configuration i. The sampling configurations are generated to obey the 
Gibbs distribution. 
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The expression of the expectation values by the weighted sampling is 
given by 

1 M 

(^sampling = ^~]A(l), (3.3) 

i=l 

where the summation runs over the sampled configurations. In the large M 
limit, (A) sampUng converges to (A). 

The configuration obeying the Gibbs distribution can be obtained from 
the Markov chain which comes from the master equation as 

where C is the transition matrix defined in the previous chapter. The matrix 
element which is the transition probability from ith state to jth state has 
a non-zero value, when the configurations of ith and jth states differ only 
by one spin. Therefore, the ith state can be modified to N states. We 
create a new configuration from the ith state by the following steps: Choose 
a site randomly and determine the direction of the spin on the site by the 
transition ratio from the configuration before flip to the configuration after 
flip. Next, repeat this procedure N times. These N trials of the spin flip 
are called "one Monte Carlo step". The Markov chain is obtained from the 
calculation of this Monte Carlo steps. In principle, every configuration can 
change to any configuration in one Monte Carlo step, because all spins can 
flip once. However, the distance between two configurations is close, if one 
configuration is obtained from the other configuration by one Monte Carlo 
step. For the precise calculation, the interval of the sampling has to be long 
enough. 

3.2 Quantum Monte Carlo Method 

In quantum systems, the definition of the expectation value is not the same 
as in the classical systems. The expectation value of a quantity A is given 
by 

TrAe-P n 
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where 7i is the Hamiltonian. To calculate the operator e~^ n is difficult for 
large-size systems because of non-diagonal terms of the Hamiltonian 7i. The 
same procedure as in the classical Monte Carlo method can not be applied 
to quantum systems, if we can not diagonalize the Hamiltonian. 

The following procedure avoids the above difficulty. Now we consider the 
the Ising system with transverse field T as a quantum system, 

ij i 

This Hamiltonian is not diagonal in the ^-representation. By the Trotter 
formula, the operator e _/3W can be described as a product of many operators 
diagonalized in z- or x-representation: 

e-W = lim ^/M e B/M )M (3J) 

M->oo 

where A = E^^j = and B = 7 of (7 = 0T). This 

decomposition is called the Suzuki- Trotter decomposition |3R |3"5|| . Using this 
decomposition, we obtain the partition function of the Mth decomposition 

Z M as 



E ({^i}l^ /M IKi})({^i}l^ /A/ l{^}) 

{o-Jfe=±l} 

x({a j2 }|e A / M |K 2 })({ ( T' }|e B ^|{a j3 }) 



x 



x({^M}|e A / M |{ ( x; A/ })({ ( x;. M }|e^|{a jl }), (3.9) 
where |{cx.,fc}) is the Mth direct product of eigenstates {cr,-} defined as 

\Wjk}} = Wji) ® 1^2) ® • • • ® I^-m)- (3.10) 

From this decomposition, the partition function is expressed as a trace of 
products of diagonalized matrices. Each part of the product is rewritten as 
follows: 

({a jk }\e A / M \{a' jk })=exp (^J^ KlJ a lk a jk \ f[6(a jk ,a> k ). (3.11) 

\ ij / 3 
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(Wjk}\e B/M \Wj,k+i}) = a-M ex P {lM ^2 ' ( 3 - 12 ) 

where 

aM = {l Shlh {ll)} ' 7M = ^ lnC ° th (^)- 
The partition function Z is represented as 

Z = lim < M Y 

Kfc=±l} 



{M / M 
[ Jf 5 KijCTikajk + lM S a 3 ka i> k +i 
fc=l \ j=l 



= lim < M V 

{<x jfc =±l} 

\ M ( N w 

exp < /3 efT ^ I JijVikVjk + r M aj k a jik+1 J > ,(3.13) 
I fc=i V ij j=i / J 

where /3 e g = /3/M and = Jm/^s- This is the partition function of a 
(d+ l)-dimensional classical spin system at the effective inverse temperature 
/? e ff. The quantum d- dimensional partition function is mapped to a (d + 1)- 
dimensional classical partition function. The representation of the classical 
(d + l)-dimension system thus can be obtained by the usual Monte Carlo 
method (see Fig. |3.1| ). 

The framework does not change, when the longitudinal field is applied. 
It produces the term —h^2 i af in the Hamiltonian. We can deal with this 
term similarly to the term — J^a^a*-. 

3.3 Path-Integral Monte Carlo 

The time development of the wave function in quantum systems is given as 

|*(r)> = Te-^o T<flW W|tf(0)), (3.14) 
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Figure 3.1: The <i-dimensional quantum system is mapped to a (d + 1)- 
dimensional classical system at the inverse temperature /3 e g. The original d- 
dimensional space is expressed as the "real space" and the additional space 
is expressed as the "Trotter". Copies of the original system without the 
quantum term are placed in the Trotter direction. The real space interactions 
are random variables in the real space direction but uniform in the Trotter 
direction. The interactions between the Trotter slices, Tm, are uniform and 
depend on the amplitude of the transverse field T. 



where the symbol T means that the operator acts on the ket vectors with time 
ordering. Decomposing the product of the short time (St = T/M) operators, 
we have 



. . . . • Tn M-l ■ Tn M-2 

\^(T)) = lim e e 

M^QO 



. TH 1 - TH 
. e l — e % M |^(0)), 



(3.15) 



where Hj = H(Tj/M). By replacing the time with the imaginary time 
(iT — > (3), we obtain 



\ty(-i8)) = lim e" 

M-»oo 



07in 

•e-^|*(0)). 



(3.16) 
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Figure 3.2: (a) A path of the imaginary time evolution is expressed as a 
configuration on a (d+ l)-dimensional space. As the transverse field becomes 
weak, the interaction between the Trotter slices, Tm(P), becomes strong. 
This fact is illustrated as the thickness of bonds, (b) The initial states |$(0)) 
are located at each end and evolute toward the center. 
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The expectation value of A is obtained as, 

(A(T)) = (A(-i(3)) (3.17) 
= (9(-ip)\A\9(-ip)) (3.18) 
= lim (*(0)\e-V M n °--- e- plM Hm ^A 

M— »oo 

xe -f3/M Hm-x . . . e -P/M WoJ^Q)). (3.19) 

If the imaginary time (3 is a real number, this formulation is similar to 
the quantum Monte Carlo method. Dividing the Hamiltonian into two parts, 
the terms including o z or a x , and inserting the complete sets between each 
product, we can map the equation ( |3.19| ) to the weighted sampling of the 
Monte Carlo simulation. One configuration in the Monte Carlo simulation 
corresponds to one possible path of the imaginary time dynamics in the 
original system. The summation of the samplings is the discrete version of 
the imaginary time path integral. The system simulated in the Monte Carlo 
is shown in Fig. [3.2| . The initial state |^(0)) is located at the right and 
the left sides. The evolved state locates at the center. A path of 

the sampling is obtained from the snapshot of the thermal equilibrium state 
of a (d + l)-dimensional space. The horizontal interaction comes from the 
transverse field. If the field is scheduled (decreases as a function of time), 
the interaction changes stronger toward the center. The same systems are 
located in the horizontal direction. These are similar to the Trotter slices in 
the quantum Monte Carlo method. 



3.4 Applying Monte Carlo Method to An- 
nealing Process 

In statistical mechanics, Monte Carlo simulation and the master equation of 
the classical system give the same result except for statistical errors. The 
results in the dynamics are also the same. (For the dynamics of the SK model, 



see [11], ^2], |43l) From this equivalence, we regard the Monte Carlo step 
as the time of the master equation. Figure |3.3| shows the probability to find 
the ground state of the SK model with N = 8 for the master equation and 
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the Monte Carlo simulation. The probability for the Monte Carlo simulation 
is obtained from the ratio of the number of the ground state configurations 
for each time in one hundred independent runs. The temperature schedule 
of both cases is T = 3/Vt. Two curves show good agreement. 

On the other hand, there is no reason that the quantum Monte Carlo 
simulation is equivalent to the solution of the Schrodinger equation. However, 
we regard that the Monte Carlo step in the quantum Monte Carlo simulation 
is also equivalent or similar to the time in the Schrodinger equation. How 
different between the two results is shown later. 

Besides the above approach, the path-integral Monte Carlo also describes 
the dynamics of the quantum system. The imaginary time path integral 
of the quantum system and the partition function at the temperature 1//3, 
where (3 is the imaginary time, have the same expression. The time develop- 
ing operator e~@ n can be decomposed like the Suzuki- Trotter decomposition 
in the quantum Monte Carlo method. Therefore, the imaginary time path 
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Figure 3.4: The results for the Schrodinger equation, the quantum Monte 
Carlo method, the imaginary time Schrodinger equation and the path- 
integral Monte Carlo method. Imaginary time Schrodinger equation means 
the results of imaginary time dynamics of the Schrodinger equation. The 
schedule of V(t) is 3/y/{t). 
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Figure 3.5: The results for the Schrodinger equation, the quantum Monte 
Carlo method, the imaginary time Schrodinger equation and the path- 
integral Monte Carlo method. The schedule of T(t) is 3/t. Except for the 
result for the Schrodinger equation, the probabilities converge to one. 
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integral is evaluated numerically. 

In this formulation, the path integral is replaced with the Monte Carlo 
simulation. However, we have three difficulties about this simulation. The 
Hamiltonian with the imaginary time 7i(—if3j/M) is no longer a Hermitian 
matrix, because the amplitude of the transverse field T(t) is a function of 
the time and generally becomes a complex number when /3 is a real number. 
We can avoid this difficulty by changing the definition of the time in T(t) 
as t ^ \ftt* , where t* is the complex of t conjugate of t, because y/tF 
is a real number for any complex number t. Secondly, the simulation is 
performed by replacing imaginary numbers with real numbers, so that we 
have to evaluate A(T) from A(J3) by the process of analytic continuation. The 
analytic continuation for the function whose values are given as a numerical 
number is a great problem itself and needs large computational power like 
the simulated annealing. Finally, the computational time of the path-integral 
Monte Carlo method may be longer than the quantum Monte Carlo method. 
The imaginary time evolution is mapped to the spatial axis, so that we 
have to wait until the system becomes in equilibrium. We can decrease the 
transverse field in the quantum Monte Carlo method without care that the 
system is trapped at a local minimum or not. In the path-integral Monte 
Carlo, we care about that the system is trapped at a local minimum or not, 
because the expectation value is obtained from the average of various paths 
whose distribution obeys the Gibbs distribution of the (d + l)-dimensional 
system. 

If we accept that the imaginary time expectation value A((3) is regarded 
as a real time one A(T), the final problem is still a difficult problem. It takes 
long time to wait until the system becomes in equilibrium, if the landscape 
of the system is complicated. 

We compared the above two approaches, the quantum and the path- 
integral Monte Carlo methods, with the results of the direct solutions of the 
Schrodinger equation. The system is the SK model with N = 8 and the 
schedule is T = 3/yi. The probability to find the ground state is plotted 



in Fig. |3.4j . The quantum Monte Carlo is performed with M = 1000. We 
wait 5000 Monte Carlo steps for the initial relaxation and take average over 
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the next 5000 steps in the path-integral Monte Carlo simulation. We put 
the inverse temperature f3 = 1000 for each calculation. The result of the 
imaginary time dynamics of the Schrodinger equation is also plotted. 

In Fig. [3.5|, we plot the solution whose schedule is V — 3/i. In this 
schedule, the solution of the Schrodinger equation does not converge to the 
ground state with probability one, but the others converge to one. The curve 
of the path-integral Monte Carlo agrees with the curve of the imaginary time 
Schrodinger equation. 

From these results, we can conclude that the quantum Monte Carlo 
and the path-integral Monte Carlo simulations are not the same as the 
Schrodinger equation. We regard that the assumptions of the two Monte 
Carlo simulations are not correct. The assumptions are the following: For 
the quantum Monte Carlo simulation, we assume that the Monte Carlo step 
equals to the real time. For the path-integral Monte Carlo simulation, we 
assume that the imaginary time dynamics equals to the real time one. 

However, the goal is to find the ground state as fast as possible. The two 
Monte Carlo methods find the ground state configuration with probability 
one, even if the schedule is too fast to solve for QA by the Schrodinger 
equation. For this purpose, we can adopt one of the Monte Carlo methods. 
We consider that the quantum Monte Carlo method has the advantage. The 
quantum Monte Carlo method is similar to SA, because the control parameter 
(T for SA and V for QA) decreases as a function of the Monte Carlo step. 
Thus, the longer we perform the simulation, the larger the probability to 
find the ground state becomes. On the other hand, a long simulation makes 
the statistical errors small for the path-integral Monte Carlo method, but 
the probability does not increase. For these reasons hereafter we adopt the 
quantum Monte Carlo method for QA of the large-size systems. 

We compared the probabilities to find the ground state in the Monte Carlo 
version of SA and QA for the SK model with N = 51 and T = T = 3/v^ 
in Fig. |3.6| . In this calculation, we perform 100 independent runs of long- 
time (t = 100000) SA with the sufficiently slow schedule T = 3/ln(t + 
1) and determine the ground state energy beforehand. Using this ground 
state energy, we can calculate the probability to find the ground state. The 
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Figure 3.6: Overlaps with the ground state are plotted. The system is the 
SK model with N = 51 and 3/y/i schedule. 

probability to find the ground state is the ratio of the number of the ground 
state configurations for each time in ten thousand independent runs of SA 
and Trotter slices of QA. The probability in QA is larger than SA. This result 
means that the QA in the quantum Monte Carlo method also improves the 
performance in finding ground state. The quantity 1 — P(t) is plotted in 
Fig. |3.7| under the log-log scale. The curve converges to zero asymptotically 
as \jt. This implies that the system almost follows the stationary state and 
the annealing schedule is sufficiently slow. Moreover, we apply fast schedule 
to the same system and the probability converges to one asymptotically as 
1/t 2 . In this fast schedule, system also follow the stationary state. 

We also plot the average energy of SA and QA under the same condi- 
tion in Fig. |3.9| . As we need the ground state of the classical term of the 
Hamiltonian, we neglect the contribution of the transverse field term in the 
calculation of the energy. We take the average energy which does not contain 
the interactions between the Trotter slices for each snapshot of the Trotter 
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Figure 3.7: The quantity 1 — P(t) is plotted. The conditions are the same 



as in Fig. [3.6 . 
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Figure 3.8: The quantity 1 — P(t) is plotted. The schedule is T = 3/t. 
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slices in QA. The average of the energy per spin is calculated by M inde- 
pendent runs for SA and M Trotter slices for QA, respectively. In spite of 
the large difference between SA and QA in P(t), the energy of QA is lower 
than SA, but the difference is not so large. This can be understood that if 
the system is trapped in a local minimum, the trapped system also lowers 
the average energy but does not count in P(t). 

Calculating the probability to find the ground state is a useful index to 
check which method has a good performance. However it is too difficult to 
perform this calculation for large-size systems, because the calculation needs 
the ground state configuration, but we do not know it. In this case, we regard 
that the method whose average energy is lower has a better performance. If 
the average energy of QA is lower than SA, QA may find the ground state 
more efficiently than SA. 

On the other hand, another index can be considered. That is the lowest 
energy in the M independent runs of SA or in M Trotter slices of QA. How- 
ever, comparing the performance by the average energy gives more precise 
conclusion than by the lowest energy, because the lowest energy depends on 
the initial state more than the average energy. If the initial state is in the 
basin of the ground state, the system goes to the ground state with large 
probability. We will calculate SA and QA by the Monte Carlo method for 
the large-size systems up to iV = 10000 and compare their performance by 
the average energy. 

3.5 Results of the Monte Carlo simulation 

Two additional parameters which do not appear in the Schrodinger equation 
are needed, when the calculation of QA is performed by the quantum Monte 
Carlo simulation. The first one is the (inverse) temperature T(= and 
second one is the Trotter number M. The original dynamics, the Schrodinger 
equation, does not contain thermal fluctuations, so that we have to take the 
limit (3 — > oo. The Trotter number M should also be infinity to take into 
account the quantum effect correctly. 

Each Trotter slice is the classical system with the effective inverse tern- 
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Figure 3.9: The average energy is plotted for SA and QA. The system is 
the SK model with N — 51 and 3/\/i schedule. 

perature f3 e g(= (3/M) except for the interactions between Trotter slices. We 
consider the parameters f3 e g and M instead of the parameters (3 and M, be- 
cause it is natural to regard that we simulate the classical (d+ l)-dimensional 
system at the effective inverse temperature. The infinity limit of the Trotter 
number also means the infinity limit of the inverse temperature, because the 
ratio (3/M is fixed. 

3.5.1 Dependence on the temperature 

First, we consider the effective inverse temperature /3 e ff- At the initial time, 
the strength of the interaction between Trotter slices is small and each slice 
is almost an independent classical system. At the end, the strength goes to 
infinity. The M spins in the Trotter direction take the same value because of 
the interactions of the infinite strength and the freedom of this direction van- 
ish. The (d+ l)-dimensional system which is equivalent to the ci-dimensional 
quantum system reduced to the <i-dimensional classical system. This can 
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be understood naturally, because the transverse field T(t) vanishes and the 
Hamiltonian goes to the classical Hamiltonian. 

The spins in the Trotter direction take same value when the interactions 
between the Trotter slices is stronger than the thermal fluctuation. If the 
temperature is high, we need long time until the interaction becomes large 
enough. From this, the temperature should be low to save the time for 
the calculation. On the other hand, the temperature is required to be high 
enough for searching the configuration space widely at the initial time. If 
the temperature is low, the system can not visit various configurations and 
trapped at a local minimum. 

To satisfy the above two conflicting conditions, we have to choose mod- 
erate values of temperature. The system can seek all possible configurations 
when the temperature is above the critical temperature of the order- disorder 
transition of the classical system. Thus, we choose /3 e e = 1 for the SK model. 

The dependence of the probability on /3 e g is shown in Fig. RTTDl The 



condition is (3 eS = 1/2, 1, 2, 5, N = 8 and M = 1000 for the SK model in the 
3 / y/i schedule. The simulation with (3 e a = 1 has a best performance among 
the four. 

3.5.2 Dependence on the number of Trotter slices 

Secondly, we consider the Trotter number M. A simulation with large num- 
ber of Trotter slices gives precise results, but the computational power for 
the calculation becomes large. The computational time for the quantum 
Monte Carlo simulation is M times longer than the conventional simulated 
annealing. We have to estimate the reasonable Trotter number of the correct 
simulation. 

We check the dependence of the simulation on the Trotter number. The 
probabilities to find the ground state versus the transverse field are plotted 



for the cases of M = 10, 20, 50, 100, 200, 500, 1000 in Fig. gJJ. The condition 
is T = 2 /y/i, N = 625 and /3 e g = 1 for the two-dimensional Edward- Anderson 
model. We find that M = 100 is large enough. 

From the above two results, we regard that the conditions of (3 e s = 1 and 
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Figure 3.10: The probability dependence on the ratio (3/M . The simulation 
with P eS = 1 has a best performance among the four. 

M = 100 give a reasonable solution. Hereafter we will use this condition. 

3.5.3 Results of large size system for the EA model 

Next, we consider systems of larger size than the previous calculations to 
demonstrate that our method can be applied to the actual problems. The 
model is the two-dimensional Edward-Anderson (EA) model. The calculation 
per Monte Carlo step for the EA model is less than the SK model. The spin 
flip is determined by the calculation of the energy gap of the spin flip. We 
have to sum up N — 1 interactions for the SK model, because all spins interact 
with each other in such an infinite range model. For the two-dimensional EA 
model on the square lattice the number of the interactions is four. Therefore, 
we can accelerate the summation in the EA model (N — l)/4 times faster 
than the SK model. The acceleration affects both of the methods SA and 
QA. 

We choose the standard condition of the calculation of this model as 
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Figure 3.11: The magnetization dependence on the Trotter number M. 
Curves of over M = 100 are large enough. 

follows; N = 625(25 2 ) with periodic boundary conditions, /3 e ff — 1, M — 100 
and the schedule T = T = c/y/i. The critical temperature T c is considered 
zero in this model. This is not the same value of the SK model. However, 
the zero temperature dynamics only seeks the configuration whose energy is 
lower than the present configuration and tends to be trapped at local minima, 
so that we use the same condition as f3 e s = 1. 

We check how the system follows the stationary (equilibrium) state which 
means equilibrium state at fixed T for SA and ground state at fixed V for 
QA. For small-size systems, we can calculate the probabilities to find the 
ground state -Psa(^) and Pqa(£) and compare with the probabilities in the 
stationary state P|^(T(t)) and Pn A (r(t)). If the two curves, the results of 
the dynamics and the statics, are the same, the system follows the stationary 
state. For large-size systems we cannot calculate the probability to find the 
ground state, because the calculation of the probability needs the ground 
state configuration beforehand and that is difficult. Nevertheless, we can 
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check whether the system follows the stationary state or not by another 
quantity. In our case, a longitudinal magnetic field is applied to the system, 
and magnetic order grows as a consequence of this field. If the system is in 
the stationary state, the magnetization is characterized by the temperature 
for SA or the transverse field for QA and does not depend on the speed of 
the schedule. The functional form of the schedule is c/ thus the speed of 
our schedule is controlled by c (slow for large c and fast for small c). The 
system can not follow the stationary state for the rapid schedules, the cases 
of small c. We check the dependence of magnetization on c. 

The magnetization versus the temperature for SA and the transverse 
field for QA is plotted in Figs. 3.12 and 3.13| . For SA, the final value of the 
magnetization (the left end of the curves) depends on c. This means that the 
schedule is too fast to follow the stationary state. This result implies that 
the system tends to be trapped at local minima in this schedule up to c = 50. 
For QA, the final value of the magnetization does not depend on c beyond 
c = 10. Thus, we can regard that the schedule Y — cj \/t with c = 10 is slow 
enough to search the ground state of the two-dimensional EA model. The 
conclusion of these two results for SA and QA is that SA does not follow the 
stationary state and QA does in the cj yt schedule. This is consistent with 
small-size analysis by differential equations (the master equation for SA and 
the Schrodinger equation for QA). 

We average out the energy for SA with M(= 100) independent runs and 
QA with M Trotter slices. The results on the average energy are plotted in 
Fig. |3.14| . Only in the middle region of the time, the energy in QA is higher 
than in SA, but finally the energy in QA goes lower than in SA. 

As shown in the previous chapter, the performance in finding the ground 
state is improved by the quantum effect. However, the quantum Monte 
Carlo uses the Suzuki- Trotter decomposition by definition. The calculation 
of the quantum Monte Carlo needs Trotter-number M times longer than the 
calculation of the classical system. Therefore, the comparison between SA 
and QA is still incomplete. The time scales of the two computations are not 
the same. One Monte Carlo step in QA takes about M times evaluations of 
spin flips than in SA. (To be precise, the calculation in QA is M{z + 2)/z 
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Figure 3.12: The dependence of magnetization on c for the conventional 
simulated annealing (SA) 



0.28 




1,11 I 


1 


c=i 


0.26 




' *'< 




C=2 - 


0.24 
0.22 






«%/<<^ 

\ % 
\ % 


C=5 
(=10 

( =20 . 

C=50 


0.2 






\ \ % 




0.18 






\ \ I 
\ \ 1 




0.16 










0.14 






\ \\\ 




0.12 






\ h 




0.1 










0.08 




.... 1 


i 





0.1 



10 



Figure 3.13: The dependence of magnetization on c for quantum annealing 
(QA) 
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Figure 3.14: Average energy for SA and QA are plotted. The schedule is 
T = T = 10/ arid the system size is N = 625. 



times longer than the calculation in SA, where z is the number of nearest 
neighbors and "2" is the number of the interactions in the Trotter direction. 
It is exactly "M times" for the infinite range model in the thermodynamic 
limit.) Taking into account this disadvantage gives a reasonable comparison 
between SA and QA. We provide a reasonable definition of the time t' of 
which the quantity is plotted and compared as a function as, 

i! = t (for SA) t' = Mt (for QA), (3.20) 

The Trotter number is 100 in our case, so that the time is t' = lOOt for QA. 



The average energy for SA and QA is plotted in Fig. |3.15| under the 
condition of such a scaling of the time. QA improves the performance in 
finding the ground state. In the calculation, we adopt that the schedule in 
SA is 100 times slower than QA, because SA could not follow the stationary 
state in the cj\ft schedule (see Fig. |3.12j ). The total Monte Carlo step of 
SA is 100 times longer than QA, because we have to use same amount of 
the computational power to compare the two methods (QA needs 100 times 
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calculations than SA). The values of the temperature and the transverse 
field at the end of the simulations are the same. We also plot the cases of 



various schedules for the temperature and the transverse field in Fig. |3.16 



and Fig. 3.171 . QA's performance is also good in the various schedules. For 
another system size, N = 10000, the average energy in QA is also lower than 
SA as seen in Fig. |3.18| . 

At the end of simulations, we cannot reach the zero temperature or zero 
transverse field. The values are quite small (~ 0.03) but not equal to zero. 
For example, the strength of the interacrion between the Trotter slices is 1.73, 
when T = 0.1/v^lO. It is not so strong compared with the effective inverse 
temperature (3 e s = 1. The thermal fluctuation still remains. To remove 
this fluctuation, we execute zero temperature dynamics after the annealing. 
Then, almost all the Trotter slices (more than 90% for the situation of the 



calculation for Fig. 3.15 ) converge to the same configuration in QA, while 
the independent runs in SA do not. 



3.5.4 Difference between the annealing and the quench- 
ing processes 

Let us take another approach to investigate the property of QA. How effi- 
ciently does the annealing process find the ground state in comparison with 
the quenching process? We calculate the average energy of the two processes, 
the annealing and the quenching, for SA and QA. The temperature and the 
transverse field are quenched at T = Y = OA/y/lO = 0.0316 ■ • ■ , and sched- 
uled from infinity to zero by a function 10/ \/t. The quenching process of the 
quantum system is almost equivalent to the case of the work by Sato et al. [B^] 
They considered the dynamics of the quenched system with r —>■ oo. The re- 
sults are shown in Fig. |3.19| . The temperature T and the transverse field T for 
annealed and quenched systems are the same at the end of the simulations. If 
the systems stay in the stationary state, the curves arrive at the same point 
at the end. The average values of energy of SA, QA and T-quenched systems 
converge to almost the same value, but the T-quenched system is trapped 
in a high-energy state. Of course the energy of the T-quenched system is 
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higher than that of QA (see Fig. |3.20| . It is an enlargement of Fig. |3.19| ). 
The behavior of the two quenched systems are quite different, though both 
of the controlled parameters T and T are quenched. The quantum system 
relaxes faster than the classical system, when the transverse field takes the 
same value as the temperature of the classical system. This is consistent 
with the result that QA reaches the stationary state in magnetization, but 
SA does not when the coefficient c changes from 2 to 50. 



3.6 Summary 

We applied the Monte Carlo method to QA. There are two ideas to translate 
the time in the Schrodinger equation to the Monte Carlo method. One is the 
path-integral Monte Carlo which regards the Trotter direction as the time 
(more precise, the imaginary time). The other is the quantum Monte Carlo 
with "the Monte Carlo step" -dependent inter- Trotter interaction. We regard 
the Monte Carlo step as the time in the Schrodinger equation. 

The path-integral Monte Carlo is performed in the imaginary time dy- 
namics. This results in the different behavior from the solution of the 
Schrodinger equation. For the quantum Monte Carlo, the result also differs 
from the result of the Schrodinger equation. Two reasons are mentioned. 
The energy dissipation occurs in the Monte Carlo simulation, while it does 
not occur in the original Schrodinger equation. Beside this, the time in the 
path-integral and the quantum Monte Carlo methods are not same as the 
time in the Schrodinger equation. However, these methods have good per- 
formance and the ground state is found, even if the schedule is too fast to 
find the ground state for the Schrodinger equation. 

We adopt the quantum Monte Carlo for large-size QA to reduce the com- 
putational time. We have to wait until the system is in equilibrium in the 
path-integral Monte Carlo and summing up sampling paths. This takes large 
computational time. 

The Trotter number we adopt is M = 100, which is large enough for 
our calculations (two-dimensional EA model with N = 625). The results are 
compared with SA. We found that QA improves a performance in finding the 
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ground state even if we take into account that the QA need M times more 
calculations than SA. Moreover, almost all the configurations of the Trotter 
slices converge to the same configuration. 

We also checked the difference between annealed and quenched systems. 
The relaxation time of the classical quenched system is quite long. On the 
other hand, the quantum quenched system does not so differ from the an- 
nealed system QA. This implies that the quantum effect accelerates relax- 
ations. Moreover, it is remarkable that the average energy of the quantum 
quenched system is lower than SA. The short relaxation time makes the de- 
pendence on the schedule small. This can be seen in Figs. [3.12 and 3.13 . 
A rapid decrease in the characteristic relaxation times for the random Ising 



magnet, LiHo .i67Yo.833F4, with the transverse field is observed jPJ]. We 
consider that the results may come from the same mechanism. 

From the application side, the optimal schedule is not known in general. 
If the dependence of the performance to find the ground state on the anneal- 
ing speed is small, we can decrease the rate not to reach the ground state 
or the optimal configuration. This is an attractive property for the actual 
applications. 
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Figure 3.15: Time evolution of the average energy for QA and SA. The 
system size is A^ = 625. One hundred Trotter slices are averaged for QA 
and 10 runs are averaged for SA. For "SA C = 2" and "QA C = 2" the 
schedules of the temperature and the transverse field are the same function. 
"SA C = 20" and "QA C = 2" take the same value of the temperature 
and the transverse field T = T = 0.1/\/T0 = 0.0316 • • • at the end of the 
simulations. 
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Figure 3.16: Time evolution of the average energy for QA and SA. The 
system size is A^ = 625. One hundred Trotter slices are averaged for QA and 
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Figure 3.17: Time evolution of the average energy for QA and SA. The 
system size is N = 625 and the schedule is T = 50/y/t and T = One 
hundred Trotter slices are averaged for QA and 10 runs are averaged for SA. 
The temperature and the transverse field at the end of the simulations are 
the same value T = T = 0.05. 
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Figure 3.18: Time evolition of the average energy for QA and SA. The 
system size is N = 10000 and the schedule is T = 20/v^ and T = 2/Vi. One 
hundred Trotter slices are averaged for QA and 10 runs are averaged for SA. 
The temperature and the transverse field at the end of the simulations are 
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Figure 3.19: Solutions for annealing and quenching processes are plotted. 
The annealing schedule is T = T = 10/ \fi and the quenched value is T = 
T = 0.1/ vTo = 0.0316 ■ ■ ■ which is the final value of the schedule. 
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Figure 3.20: Solutions for annealing and quench process for QA are plotted. 
This is the enlargement of Fig. |3.19. 



Chapter 4 

Application to the Traveling 
Salesman Problem 



All the problems we calculated up to the previous chapter were on the Ising- 
spin systems. However, the combinatorial optimization problem is not lim- 
ited to the problems based on physics. One of the major problems is the 
traveling salesman problem (TSP) [|I], |2[. We will apply our method to the 
TSP in this chapter. The TSP is a hard problem to solve by a computer, 
say an iVP-hard problem. Therefore, comparison between SA and QA in the 
TSP may provide a good measure of their performances in the problems of 
other areas not limited to physics. In other words, we will check generality 
of QA for combinatorial problems. 



4.1 The Traveling Salesman Problem 

We have N points or cities in a space with distances dij between them. The 
task is to find the minimum-length closed tour that visits each city once 



and returns to its starting point. Figure |4.1| shows an example. Practical 
applications include scheduling of truck deliveries, airline crews, and the 
movements of an automatic drill-press or robot arm. In many applications 
the "distances" dy are abstract quantities not related to a Euclidean dis- 
tance between points in a real space; they may not even satisfy the triangle 
inequality d^ < d^ + dy. For the example of truck deliveries, the cost to 
transport freight from city i to j is affected by many factors, amount of the 
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Figure 4.1: The traveling salesman problem, showing a good (a) and a bad 
(b) solution to the same problem. 

freight to city j, labor, gas, highway charge, the weather, etc. The length of 
the tour is defined by 



N 



5> 

m=l 



i(m),i(m+l) j 



(4.1) 



where i(m) means the city of the mth stop, the tour is closed as i(N + 1) = 
and N is the number of cities. 

One can think of a simple approach, listing all possible tours of the prob- 
lem to search the optimal tour. The number of all the possible tours is 
(N — l)!/2. This number increases rapidly and becomes over the limit of the 
computer power. Thus, a faster algorithm is needed to solve the actual prob- 
lems. Some algorithms are optimized for TSP Jl|, Q and others are general 



algorithms which can be applied to TSP pl| . The simulated annealing and 
the quantum annealing belong to the latter. 

4.2 Simulated Annealing on TSP 

The simulated annealing is based on the Monte Carlo simulation, so that we 
have to determine the dynamics, how to rearrange of the tour. We choose 
the simplest rearrangement of the tour like one-spin flip in the Ising-spin 
system. That is to exchange two stops of cities i and j, which is illustrated 
in Fig. O. 
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This rearrangement of the tour is realized in the following way: (1) Choose 
two cities i and j. (2) Exchange stops of cities i and j with the probability 



wa->b 



exp(-/3L B ) 



exp(-(3L A ) + exp(-(3L B ) 



(4.2) 



where A is the route before exchange and B is the route after exchange. 

One Monte Carlo step consists of N(N — l)/2 exchange trials of all pairs. 
We fix the starting point and reduce the trials to (N — 1)(N — 2)/2, because 
the tour is closed. The reversibility of the tour still remains in our calculation. 
The temperature T is scheduled in a form cj y/i, where t is the Monte Carlo 
step. 





J 



J 



Figure 4.2: The simplest rearrangement of the tour. This rearrangement is 
expressed in four spin-flip dynamics. 



4.3 Quantum Annealing on TSP 

We have to introduce a quantum fluctuation to the TSP, but that is not 
trivial. The system is not the physical system. We cannot easily imagine 
how to introduce a quantum fluctuation. A straightforward way is that we 
express TSP as a binary system namely, the Ising system, so that we can 
introduce the quantum effect. The idea to express TSP as a binary system 



was proposed by Hopfield and Tank H5L E6 . 
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We choose stochastic binary units n ia to represent possible solutions: 
rii a = 1 if and only if city i is the ath stop on the tour and n ia = for 
otherwise. There are N 2 units altogether. The spin representation can be 
obtained by the transformation rii a — > (<7j a + l)/2. The total length of the 
tour is 



2 

ij,a 
ij,a 

and there are two constraints: 



L = - ^ dijUjajn^g+i + % ;0 -i) (4.3) 
njtfiofea+i + ff j,a-i) + const. , (4.4) 



^^riia = 1 (for every city i), (4.5) 

a 

^^n ia = 1 (for every stop a). (4.6) 

i 

The first constraint says that each city appears only once on the tour, while 
the second says that each stop on the tour is at just one city. The tour is 



expressed as a configuration of N x N units as illustrated in Fig. [O 

The simplest rearrangement of the route of the traveling salesman is ex- 
changing the stop of cities i and j. In the spin expression, four spins are 
flipped at a time. The salesman stops at cities i and j for the ath and the 
6th stops respectively. On- route spins take a ia = = 1, and off-route spins 
take (Tib = <Jja — — 1- One up spin and one down spin in the same row are 
flipped at a time, so that the summation which runs in the horizontal direc- 
tion does not change. In the vertical direction, the same situation occurs. 
The constraints ( |4.5|) and ( |4.6|) are kept in the four-spin flip. The Monte 
Carlo simulation in the spin expression is performed by this procedure. 

We fix an — 1 which means that the tour starts from city 1. Possible ex- 
changes with any two stops of cities are the combination number of choosing 
two from iV — 1 as (N — 1)(N — 2)/2. One Monte Carlo step must contain 
(N — 1)(N — 2)/2 spin-flip trials. This number of trials is suitable that the 
system has iV 2 spins and 0(N 2 ) trials are performed in one Monte Carlo 
step like the usual one-spin flip procedure. 
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Figure 4.3: The illustrated tour is3->2-^4^1 with iV = 4 TSP. 
Solid and open circles denote units that are "1" and "0" respectively. The 
condition is open boundary. 



Let us consider to introduce quantum effects to the classical system. The 
dynamics of SA was described by one spin-flip procedure in the preceding 
chapters, and we chose the quantum effect which flips one spin at a time. 
That is the transverse field term, — T^2 i af. The linear combination of the 
off-diagonal operator of has a non-zero value, when a bra state differ from 
a ket state only by one spin. If the two spin-flip dynamics is applied to SA, 
we have to choose the square term of erf as — T £\ . erf a*. As known in the 
Suzuki- Trotter decomposition of the XY or Heisenberg models, the two-spin 
interaction is mapped to the four-spin interaction in the classical Ising model 
of the (d + l)-dimensional system |47] . 



On the other hand, the traveling salesman problem mapped to the Ising 
system has to flip four spins at a time. The quantum fluctuation term should 
be the four-spin interaction term as — r Yliijki a i cr j (J k 'i '■ We have an eight- 
spin interaction term by the Suzuki- Trotter decomposition of the Hamilto- 
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nian with this term. Figure [4.4] shows the interactions in a decomposed 
classical system of one-, two- and four-spin interactions of the quantum sys- 
tem. Spins which contact through the gray domain in the figure interact 
with each other. The computational cost is expected to increase when the 
number of spins which interact each other increases. Moreover, various types 
of spin-flip dynamics have to be introduced to recover ergodicity for such a 
system. 



: : : : 



: 




j 



j 



k 



Suzuki-Trotter decomposition 



(7) (i j) Q j k T) 

Figure 4.4: By the Suzuki- Trotter decomposition one-, two- and four-spin 
interactions transform to interactions which contain spins twice as many - 
two, four and eight spins. 

To avoid such difficulties, we adopt the same quantum fluctuation we 
used before as — T af. The same approach as in the previous chapter can 
be applied. The system is mapped to a (d + l)-dimensional system by the 
Suzuki- Trotter decomposition. However, the spin-flip procedure of QA for 
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TSP differs from the procedure of QA in the previous chapter. The four 
spins in a Trotter slice are flipped at a time to keep the constraints, though 
one spin was flipped so far. The spins are flipped by obeying the transition 
probability which is obtained from the energy gap between the states before 
flip and after flip. We computed the energy gap from the classical energy 
term of the Hamiltonian (the tour length) and the interactions between the 



Trotter slices for the four spins. That is shown in Fig. [4.5| . We can perform 
QA for the traveling salesman problem in this procedure. 




Figure 4.5: The four spins (big circles) are flipped at a time. The transition 
probability is computed from the length of the tour and the two-spin inter- 
actions of the nearest neighbors in the Trotter directions, (the interactions 
between big and small circles) 



4.4 Results of Monte Carlo Simulations 



We perform SA and QA for the four cases of TSP whose optimal tours are 
shown in Fig. [4.6| . The first problem "random" is the case in which the cities 
are located at random. Cities are located at random but clumped with dense 
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and spares regions for the second problem "semi-random" . The optimal route 
of the third problem has an "H" shape, (we call this problem "H-character" 



hereafter.) The fourth problem is "ulyssesl6" of TSPLIB P8| . The number 
of cities is N = 16 for all the problems. 

The cities are located on a square with the side length to make the 
length of the tour extensive for "random" , "semi-random" and "H-character" . 
For "ulyssesl6", we re-scale d^ and set the average to 2.2. The average, 
the dispersion and the ratio of the dispersion and the average are shown in 
Table O 




Figure 4.6: Optimal tours of the four problems. 



Simulations are performed with 500 Trotter slices for QA and 100 inde- 
pendent runs for SA. We observe two quantities, the probability to find the 
minimum-length P(t) and the average of length (L). The probability P(t) is 
obtained from the ratio of the number of the ground state configurations for 
each time in 100 independent runs for SA and 500 Trotter slice for QA. The 
probabilities of SA and QA for the four problems with the 10/ \/i scheduling 
are plotted in Fig. |4.7| and the energy values are in Fig. |4.8| . The plots are 
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random 


semi-random 


H-character 


ulyssesl6 


average (a) 


2.195 


2.642 


2.155 


2.200 


dispersion (b) 


0.973 


1.152 


0.867 


1.556 


b/a 


0.443 


0.436 


0.402 


0.707 



Table 4.1: The average, the dispersion and the ratio of the dispersion and 
the average for the four problems. 



located in the same order as Fig. [4.6| . In this schedule, QA has a better 
performance than SA for all the four problems. Moreover, we also calculated 
for the case of 5/v^ and 2/y/i schedules (see Fig. |OH4.12| ). Results show 
that QA is better than SA in both of the average length and the probability 
for all problems. The probabilities of SA for the "H-character" problem are 
obviously saturated and fall short of probability one, and the probability of 
QA almost reaches one. This is similar to the situation in the differential 
equations for small-size systems for SA and QA. 

Next, we investigate the asymptotic property of finding the optimal state. 
The quantity 1—P(t) is plotted in Fig. |4.13| under the log-log scale. The curve 
is proportional to 1/t in the asymptotic region. This implies that the system 
almost follows the stationary state and the annealing schedule is sufficiently 
slow. The reason is discussed in Sec. |2j]. 

Lastly, we compare the annealed and quenched systems of SA and QA 
to make the effect of the annealing clear. Annealing is introduced to avoid 
slow relaxation in low or zero temperature in the system whose landscape is 
complicated. If the annealing works, the length for the annealed system is 
lower than the quenched system. The quenching parameter is the tempera- 
ture T for SA and the transverse field T for QA. These values are quenched 
to 0.1/vT0 = 0.0316 ■ ■ ■ . In the annealed system, the temperature and the 
transverse field decrease from infinity to the same value as in the quenched 
system. We calculate for the two problems, "random" and "ulyssesl6" . The 
results of the probability and the average energy are shown in Fig. [4.14 



and |4.15| . We can list the four system in the order of the performance to 
find the ground state: (1) the quantum system with annealing "QA", (2) 
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Figure 4.8: Average of the length. The scheduling is 10/y/i. 
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Figure 4.9: Probability to find the minimum-length. The scheduling is 
5/Vt. 
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Figure 4.10: Average of the length. The scheduling is 5/y/i. 
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Figure 4.12: Average of the length. The scheduling is 2/y/i. 
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Figure 4.13: The quantity 1 — P(t) is plotted. The transverse field is 
scheduled as T = 10/y/i. 

the quantum system without annealing "QA(quench)", (3) the classical sys- 
tem with annealing "SA" and (4) the classical system without annealing 
"SA(quench)". We can find that both of the annealed systems show better 
performance than the quenched systems. It is notable that the probability 
for the quenched quantum system "QA(quench)" is larger than SA. In spite 
of annealing the temperature, SA finds the ground state less often than the 
quenched quantum system. From these results, we obtain the following con- 
clusions. The relaxation of the quantum systems is faster than the classical 
systems whatever we control the parameters, annealing or quenching. The 
annealing process accelerate the relaxations for both of the quantum and the 
classical systems. 
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4.5 Summary 

We have applied QA to TSP by mapping the problem to the Ising system. 
SA has a four-spin flip dynamics in the spin representation of the TSP. QA 
should contain the quantum effect which flips four spins at a time. Because 
of the difficulty of the calculation, we have adopted the one-spin dynamics 
by applying the transverse field. 

Numerical results show that QA has better performance than SA in the 
probability to find the minimum-length of the closed tour and the average 
of the tour length and obeys the 1/t-law asymptotically. We only calculated 
the case in which the Trotter number is 500 for QA. The dependence on the 
Trotter number was not studied, and we did not estimate how many Trotter 
slices the quantum Monte Carlo simulation need. We also did not check 
the actual performance of QA. QA needs the computational power Trotter- 
number times as large as the calculation of SA, but we did not take into 
account this disadvantage. We only calculated the quantities as a function 
of the Monte Carlo steps, not the real computational time. In other words, 
we cannot conclude which methods find the minimum-length tour of TSP 
in small computer power, SA or QA, and only conclude so far that QA 
can be applied to TSP and QA improves in terms of the Monte Carlo step 
not the real computational time. These are the future problems. Another 
future problem is to study the effect of the four-spin interaction term as a 
quantum effect. It is not clear that the multi-spin interaction makes the 
system converge to the optimal state faster or not. 

If QA for TSP accelerates the search for the optimal solution actually, 
QA may be applied to various problems which can be mapped to the Ising- 
spin systems. We should have used a four-spin transition term as a quantum 
effect in QA for TSP to have correspondence with the four-spin flip dynamics 
of SA, but we adopt a one-spin transition term — T £\ of. The results show 
that QA finds the optimal solution and the one-spin transition term works. 
This term forces neighboring Trotter slices to be in the same configuration 
in the representation of Suzuki- Trotter decomposition. Assuming that this 
effect is the essence of QA, we can import such an effect to non-Ising systems. 
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For instance, it is possible to extend our method to the problems represented 
in terms of Potts spins. 
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Figure 4.14: Difference in the probability between annealed and quenched 
systems. Top and bottom figures are the results of "random" and "ulyssesl6" 
respectively. 
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Figure 4.15: Difference in the average energy between annealed and 
quenched systems. Top and bottom figures are the results of "random" and 
"ulysses!6" respectively. 



Chapter 5 
Summary 



We have proposed a new idea on the general method to solve combinato- 
rial optimization problems. The idea is quantum annealing (QA) in which 
quantum tunneling effects cause transitions between states, in contrast to the 
usual thermal transitions in simulated annealing (SA). First, we studied QA 
by the time-dependent Schrodinger equation for small-size systems. Next, we 
extended QA to the stochastic process, the quantum Monte Carlo method. 
The two dynamics, the Schrodinger equation and the quantum Monte Carlo 
dynamics, are not equivalent to each other. However, for the purpose of find- 
ing the ground state, both of them work and the performance is improved 
in comparison with SA. Finally, we presented the possibility of the imple- 
mentation of QA for the general combinatorial optimization problems. The 
test-bed was the traveling salesman problem (TSP). The chance of finding 
the optimal solution by QA was larger than by SA. 

The idea of QA is based on SA which was proposed by Kirkpatrick et 
al. to find the optimal state of the optimization problems by introducing the 
annealing process [|2T|. Annealing is performed by decreasing the temperature 
of the system from a high temperature to a low temperature. For instance, 
let us consider the silicon oxide Si02- If the decreasing rate is slow enough, 
the system is in crystalline state. If the rate is fast, the glass state appears. 
SA imitates this process by the Monte Carlo simulation. We consider QA by 
replacing thermal fluctuations in SA with quantum fluctuations. 

First, we adopt the transverse field as the quantum effect of the Ising 
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spin systems. The task is to find the ground state of the system without 
the transverse field. The transverse field term is controlled and vanishes 
in the infinite-lime limit so that the system approaches the ground state. 
The dynamics of this system is not described by the Monte Carlo simulation 
used for SA. The natural dynamics of the present system is provided by 
the Schrodinger equation. We checked the performance of QA in finding 
the ground state in comparison with SA by numerical calculations. In these 
calculations, we applied the same function to the time- dependent schedules 
of the temperature T and the amplitude of the transverse field T. For the 
cj \ft schedule, the behaviors of QA and SA were quite different. By QA we 
find the ground state with probability one, while we do not find the ground 
state by SA. We also found the remarkable property that the probability to 
find the ground state converges to one asymptotically as 1/t. This property 
suggests that the system follows the stationary (or almost stationary) state, 
because the probability to find the system in the ground state is expressed 
using the perturbation in terms of 1) as 



where £7^ is the energy of the ith state of the non-perturbed (classical) 
system and E^ is the ground-state energy. When we put T — c/ \Jt, the 
1/t-law appears. 

Secondly, we have considered the calculation for the large-size systems. 
The calculation of the Schrodinger equation faces difficulties for larger N 
because the number of states increases exponentially as 2^ and the size of 
the storage for the Hamiltonian reaches the storage limit of the computer. 
The Monte Carlo method can be used to avoid this difficulty There are two 
ideas to replace the dynamics of the Schrodinger equation with the Monte 
Carlo method. One is the path-integral Monte Carlo Q which regards 
the Trotter direction as the time (more preciously, the imaginary time). The 
other is the quantum Monte Carlo with "the Monte Carlo step" -dependent 
interaction between Trotter slices. We regard the Monte Carlo step as the 
time in the Schrodinger equation. From the similarity to SA, we adopted the 
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quantum Monte Carlo for large-size systems of QA. 

We found that QA improves a performance in finding the ground state 
even if we take into account that QA needs M times more calculations than 
SA. The Trotter number we adopt is M = 100, which is large enough for 
our calculations (for the two-dimensional EA model with N = 625). The 
1/t convergence of the probability to find the ground state still appears in 
the Monte Carlo calculation of QA. We also found that the relaxation of 
the quantum system is faster than the classical system, if the amplitude of 
the transverse field T in the quantum system and the temperature in the 
classical system have the same value. This implies that the quantum effect 
accelerates relaxations. We consider that this fact is the reason for the better 
performance of QA than SA. 

Lastly, we have applied QA to TSP by expressing the problem in terms of 
the Ising system. Numerical results show that QA has a better performance 
than SA in the probability to find the minimum-length of the closed tour and 
the average of the tour length at the same Monte Carlo step. Moreover, the 
1/t convergence of the probability to find the optimal tour still appears in 
TSP by QA. These results imply that the QA works also in TSP. However, 
we did not check the actual performance of QA. QA needs the computational 
power Trotter-number times in comparison with SA, but we did not take into 
account this disadvantage. We only calculated the quantities as a function 
of the Monte Carlo step, not the real computational time. In other words, 
we can not conclude which method finds the minimum-length tour of TSP 
with smaller computer power, SA or QA, and only conclude so far that QA 
can be applied to TSP and QA improves in terms of Monte Carlo steps, not 
in the real computational time. 

The findings in our study are summarized briefly as follows: QA provides 
a general and faster algorithm to find the ground state of the Ising system 
whose energy landscape is complicated. If we adopt the schedule of the 
transverse field as T = c/ y/i, the asymptotic form of the probability to miss 
the ground state is proportional to 1/t. The relaxation of the quantum 
system is faster than the classical system under the condition of the control 
parameter T = T. From the application side, the optimal schedule is not 
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known in general. If the relaxation is fast, we do not have to take into 
account the possibility to be stuck in local minima more seriously in QA 
than SA. QA may also work for the combinatorial optimization problems 
which can be mapped to the Ising spin system. 

The future problems are the followings: The analytical foundations of 
the better performance of QA than SA and the critical rate of the decreasing 
the transverse field are needed. (The same method of the analysis of SA 
by Geman and Geman |22]| may be applied to the analysis of the Trotter 
decomposed system of QA.) The application to the study of the ground 
state property of the Ising spin systems, especially the three-dimensional EA 
model, is attractive. Whether the same picture as the SK model, the Parisi 
picture [13, HJ, can be applied to the three-dimensional EA model or not is 
an open problem. For some problems in SA the multi-spin flip is performed 
and the corresponding interaction term representing quantum effect in QA 
is the multi-spin interaction term. For instance, SA for TSP needs four- 
spin flip, but we adopt only a single-spin interaction term (the transverse 
field). To study the effect of the multi-spin interaction term is needed. It is 
not clear whether the multi-spin interaction makes the system converge to 
the optimal state faster or not. On the other hand, the transverse field is 
mapped to the interaction between the Trotter slices and the strength of the 
interaction is controlled by the strength of the field. This interaction forces 
neighbor Trotter slices to be in the same configuration in the representation 
of Suzuki-Trotter decomposition. Assuming that this effect is the essence of 
QA, we can import such an effect to a wide class of non-Ising systems. For 
example, it is possible to extend our method to the problems represented 
in terms of Potts spins. This is one possible way to apply QA to general 
problems which can not be expressed in terms of the Ising spin. 



Appendix A 
Single-spin problem 



In this Appendix we explain some technical aspects to derive the exact so- 
lution of the time- dependent Schrodinger equation for the transverse Ising 
model with a single spin. The three cases of T(t) = —ct,c/t and c/yi will 
be discussed. 



A.l Case of T(t) = — ct (Landau-Zener model) 

m, mm, mm 

Let us express the solution of the Schrodinger equation at time t by the 
parameters a = (+\if>(t)) and b = ( — \ip(t)). The Schrodinger equation Q2.3Q 
with 7i = —ha z — Ta x is expressed as a set of first order differential equations 
for a and b. It is convenient to change the variables as 

a=^=(a + b), b = -=(a-b), (A.l) 

by which the Schrodinger equation is now 



d 2 b(t) 
dt 2 

By using the notation 



+ (-zc + h 2 + c 2 t 2 )b(t) = 0. (A.2) 



z = V2^e- m/ H, (A.3) 
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we find 

d 2 b(t) 1 1 . , , . 

+ (A + ---**)&(*) = 0. (A.5) 

The initial state is specified as a = b = l/y/2 or b = as t — ► — oo. The 
solution of (|A.5 ) satisfying this condition is the parabolic cylinder function 



iz) |5T|. Thus, we obtain the solution as 



m = i^-^my (A . 6) 

b(t) = CtD-x-idz), (A.7) 
where C\ is a constant. To fix C±, we use the condition 

|5(-oo)|= 1 =1. (A.8) 

Then we have 

d = -5Le~^ 2 / 8c . (A.9) 
v 7 ^ 



The wave function of this system is given in Eq. ( |2.6|) . 

A.2 Case of Y{t) = c/t 

We next consider the case of T(t) = c/t. By eliminating a from the Schrodinger 
equation, we obtain 

d 2 b(t) 1 dT(t)db{t) 
~dt 2 T(t) dt dT 

+ { h2 ^-m^) m ^ (A ' 10) 

Substituting F(t) = c/t, we have 

d 2 b(t) ldb(t) ift c 2 \ I/N 

^ +T + ^ 6 (i) =0- (A.ll) 



cit 2 i dt V t t 2 
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The solutions of this equation are expressed by the confluent P function 

oo 

P \ ih 1 ic t 
—ih — ic 




(A.12) 



the right-hand side of which has two independent expressions in terms of the 
confluent hypergeometric function 

f(t) = e iht t ic F{l + ic,l + 2ic;-2iht), (A.13) 
g (t) = e iht t ic {-2iht)- 2lc 

xF(l-ic,l-2ic;-2iht). (A.14) 

The general solution is b(t) = C\f(t) + C2g(t). Using the initial condition 

b(0) = Cl + C 2 =-L, (A.15) 

a(0) = C x -C 2 =±=, (A.16) 

we find 

b{t) = --=e iht t ic F(l + ic,l + 2ic; -2iht). (A.17) 
v2 

The asymptotic forms of b(t) and |&(t)| 2 are then given as 

sf2{2h)- ic T{2ic) 



b(t) 



r(ic) 

x { e ' iht -^/ 2 + ce iht+ * c/2 (2ht)- 1 } , (A. 18) 



2 sinh(Trc) f ccos(2ht) cV c \ 

im ^^M^c){ + —fa— + 4Jwr (A ' 19) 
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A.3 Case of V(t) = c/y/i 



The final solvable model has T(t) = c/yi. The Schrodinger equation ( |A.10|) 
is then expressed as 

d 2 bit) ldb(t) 2c 2 + 2/A 7/N , k . 

W + (h 2 + b{t) = 0. (A.20) 



dt 2 2t dt V 2t 



The solution is the confluent P function []5l 

oo 



P 



zfo | — 27 i 



1 
2 

00 



e«*P\ -2^t >, (A.2i; 



1 \ 7 I 

where 7 = c 2 /2h. The two independent solutions are thus []5T| 

/(*) = e^F(i-i 7 ^;-2iM), (A.22) 

g (t) = e m (-2iht) 1 ^F(l-i 1 ,-;-2iht). (A.23) 

The general solution of ( |A.20|) is therefore the linear combination of the above 
two functions 

b(t) = C 1 f(t)+C 2 g(t). (A.24) 
The constants C\ and C2 are fixed by the requirement 

6(0) = = (A.25) 

O (0) = ^-ft = -L. (A.26) 
Substituting G\ and C 2 into Eq. ( |A.24j ), we find 
bit) = -Le iht F Q - vy, I -2iht 
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The asymptotic form is 



b{t) ~ \/ne 



-irc 2 /4h 



e- M (2hty 



-vy 



V2r(i-z 7 ) ' cr(-z 7 ) 



+ 



+e lhi (2/it) 



-1/2+17 



-(l/4)iri 



{ \/2r(i 7 ) ^ 2v^r(| + z 7 ) 



(A.28) 



The probability 1 6(oo) | 2 that the system remains in the excited state can be 
calculated as the asymptotic form of ( [A.28Q with the condition c 2 //i > 1 

2 



|6(oo) 



"77T 



| 7 -l/2 e (5/4)7ri 



r(|-i 7 ) r(-i 7 ) 



"77T 



"77T 



,1/2 



(-17) 



*7 



_ Z7 ] +7 -l/2 e (5/4W ( _ Z7)J7 +l/2 



'7 



47 



/l 2 



256 7 64c 4 ' 



(A.29) 

2 

(A.30) 
(A.31) 



A. 4 Dependence of the final value on the ini- 
tial condition 



We show that we can choose the initial condition so that the final state is 
the ground state when Y = c/yt. This is not possible for T = —ct or c/t. 
From Eq. ( |A.24| ), the asymptotic form of the solution as t — > oo is 



b(t) ~ C x 



^ e -^/4h-iht(2ht)- ic2 / 



2h 



T(l/2 - ic 2 /2h) 



-irc 2 /4h—iht 



(2ht) 



-ic 2 /2h 



c 2 T(-ic 2 /2h) 



(A.32) 



The coefficients C\ and C<i are fixed under the conditions fo(oo) = and 



+ |6(0)| 2 = 1 as 
Ci 

c 2 



smh(7Tc 2 /h) 
+ 2sinh 2 (vrc 2 /2/i) 
ic 2 T{-ic 2 /2h) 
hT{l/2-ic 2 /2h) l ' 



-1/2 



(A.33) 
(A.34) 
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This solution is not the ground state of the Hamiltonian 7Y(0). 

The reason why one cannot obtain such a solution for the other schedules 
(T = c/t, —ct) is the following: The general solution for T = c/t also has two 
coefficients, and the initial state is represented as the linear combination of 
two terms whose phases are indefinite: 

a(0) = C^X^-C^X^, (A.35) 
6(0) = C x t ic \ t ^+ C 2 r ic \ t ^. (A.36) 

The lowest-energy state at t — corresponds to a(0) = 6(0) = 1/V2 (times 
an arbitrary phase factor), which is realized by choosing C<i = in Eqs. 
( |A.35| ) and ( |A.36| ). The indefiniteness of t %c as t — ► is irrelevant because 
this is only the overall phase factor. Such a situation does not happen for 
other values of a(0) and 6(0), leading to a serious difficulty to determine the 
wave function at t — 0. Thus we cannot choose an initial condition other 
than a(0) = 6(0) = 1/V2. A similar fact exists in the case of T = —ct. 
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